Draft version August 30, 2011 

Preprint typeset using L^-T^X style cmulateapj v. 11/10/09 



GYROSCOPIC PUMPING IN THE SOLAR NEAR-SURFACE SHEAR LAYER 



o 
< 

(N 
P? 

6 



(N 
> 

o 
o 



Mark S. Miesch 1 ' 2 and Bradley W. Hindman 2 

1 High Altitude Observatory, NCAR*, Boulder, CO, 80307-3000, USA: miesch@ucar.edu and 
2 JILA and Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO, 80309-0440, USA 

Draft version August 30, 2011 

ABSTRACT 

We use global and local helioseismic inversions to explore the prevailing dynamical balances in the 
solar Near-Surface Shear Layer (NSSL). The differential rotation and meridional circulation are inti- 
mately linked, with a common origin in the turbulent stresses of the upper solar convection zone. The 
existence and structure of the NSSL cannot be attributed solely to the conservation of angular mo- 
mentum by solar surface convection, as is often supposed. Rather, the turbulent angular momentum 
transport accounts for the poleward meridional flow while the often overlooked meridional force bal- 
ance is required to maintain the mid-latitude rotational shear. We suggest that the base of the NSSL 
is marked by a transition from baroclinic to turbulent stresses in the meridional plane which suppress 
Coriolis-induced circulations that would otherwise establish a cylindrical rotation profile. The turbu- 
lent angular momentum transport must be non-diffusive and directed radially inward. Inferred mean 
flows are consistent with the idea that turbulent convection tends to mix angular momentum but 
only if the mixing efficiency is inhomogeneous and/or anisotropic. The latitudinal and longitudinal 
components of the estimated turbulent transport are comparable in amplitude and about an order 
of magnitude larger than the vertical component. We estimate that it requires 2-4% of the solar 
luminosity to maintain the solar NSSL against the inertia of the mean flow. Most of this energy is 
associated with the turbulent transport of angular momentum out of the layer, with a spin-down time 
scale of r-j 600 days. We also address implications of these results for numerical modeling of the NSSL. 



1. INTRODUCTION 

Helioseismic inversions of global acoustic oscillation 
frequencies reveal two striking rotational boundary lay- 
ers at the upper and lower edge s of the solar convective 
envelope ^Thompson et all 120031: lHowell2009D . Through- 
out the bulk of the convection zone, the angular velocity 
fi decreases by about 30% between the equator and lati- 
tudes of ±75°, with conical isosurfaces such that the gra- 
dient is primarily latitudinal. Substantial radial gradi- 
ents in n occur only near the base of the convection zone 
(0.69i? < r < 0.72i?, with R as the solar radius) where 
the super-adiabatic stratification of the envelope meets 
the sub-adiabatic stratification of the radiative interior, 
and in the surface layers (0.95i? < r < R) where deep 
convection meets the heirarchy of smaller-scale convec- 
tive mo tions sustained by rad iative cooling in the photo- 
sphere ([Nordlund et al.|[2009l ). The lower boundary layer 
is known as the solar tachocline and is thought to play 
an essential role in regulating the dynamical coupling be- 
tween the convection zone and the radiative interior and 
in generating the large-s cale magnetic fields that underlie 
the solar activity cycle (|Hughes et al.l 120071 ) . The upper 
boundary layer, known as the Near-Surface Shear Layer 
(NSSL), is similarly complex and enigmatic. Its exis- 
tence must arise from nonlinear feedbacks among turbu- 
lent convective motions spanning vastly disparate spa- 
tial and temporal scales with correspondingly different 
sensitivities to the rotation, density stratification, and 
spherical geometry. 

Although somewhat less celebrated than the 
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tachocline, the NSSL is more accessible to helio- 
seismic probing and as such, provides a unique window 
into the dynamics of the solar convection zone. The 
higher resolution of the helioseismic inversion kernals 
near the surface allows for a more reliable determination 
of the rotational gradients. Even more significantly, local 
helioseismic inversions enable a determination of the 
meridional flow throughout much of the NSSL to a de- 
gree that is not possible in the tachocline. Together with 
photospheric Doppler and tracer measurements, such in- 
versions indicate a meridional flow that is highly variable 
but systematica l ly poleward (iSnodgrass fc Dailevl 119961: 
Hathawavl 119961: lHaber et al.l 120021: iZhao fc Kosovichevl 
200l" iGonzalez-Hernandez et al 
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Hathawavl (|2011f l has recently suggested that the dy- 
namics that give rise to the poleward meridional flow 
at the solar surface occur entirely within the NSSL. This 
conclusion is based on estimates for the subsurface merid- 
ional flow obtained from correlation tracking of surface 
features. These suggest that the poleward surface flow 
reverses near the base of the NSSL, with equatorward 
counter-flow at a depth of 35 Mm below the photosphere 
(r m 0.95R). This is in contrast to previous estimates 
based on local helioseismology, which suggest that the 
poleward flow persists throughout the upper convection 
zone, with the return flow requir ed by mass conservata- 
tion o ccuring well below 0.95 R (Gil es et alJll997t iGiles] 
1999L iChou fc D ai 2001; iBeck et al.l 120021 : iBraun fc Fan! 
If confirmed, Hathaway attributes this relatively 
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shallow meridional flow structure to supergranulation. 

More specifically, the existence of the NSSL has 
long been attributed to the tendency for photo- 
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Fig. 1. — (a) Rotation rate H and (ft) specific angular momen- 
tum C inferred from global hclioscismic inversions. Note the sharp 
decrease in Q near the surface which defines the NSSL. These re- 
sults are based on RLS inversions of GONG data fro m four non- 
overlapping intervals in 1996, provided by R. Howe IIHowe et al.l 
2000; Scho u et al.1120021) . 

spheric convection (on scales of supergranulation and 
smaller) to conserve angular momentum locally, with 
fluid parcels spinning up and spinning down as they 
move t oward and away from t he rotation axis respec- 
tively (IFoukal fc Jokipfil Il975l: IGilman fc Foukall [19791: 



] thawavl |1982| iDeRosa et al.l 12002b lAugustson et al l 
: lHathawavll2011[ ).~ Such behavior is often found in 
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numerical simulations of rotating convection when the 
rotational influence is weak, such that the convective 
turnover time s cale is much less than the rotation period 
dGilmanl[l977|: IGilman fc Foukall [1979^ IHathawavl 119821: 

Augustson et al.1 
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layers where turnover time scales associated with gran- 
ulation and supergranulation are of order a day or less 
and the average rotation period is 28 days. 

Although the conservation of angular momentum is 
a valid interpretation of the numerical experiments, it 
alone cannot account for the existence of the NSSL; there 
must be more to the story. As we will demonstrate here, 
meridional forces play an essential role in determining 
the angular velocity profile in the NSSL, regardless of 
the nature of the convective angular momentum trans- 
port. Furthermore, as we will also demonstrate, the an- 
gular velocity profile of the NSSL inferred from helioseis- 
mic inversions is not consistent with angular momentum 
homcgenization. 

In this paper we propose that the poleward meridional 
flow in the NSSL is closely linked to the inward f2 gra- 
dient. The physical mechanism underlying this link is 
what we refer to as gyroscopic pumping, whereby a zonal 
forcing (axial torque) induces a meridional flow as a con- 
sequence of dynamical equilibration mediated by the in- 
ertia of the mean flows (i.e. the Coriolis force). In $2] 
and S}3]we discuss the dynamical balances that are likely 
to be achieved in the NSSL and in <g] we exploit these 
dynamical balances in order to estimate the characteris- 
tics and efficiency of turbulent transport. Then, in $5]we 
consider these results from the perspective of two simple 
theoretical paradigms. We summarize our results and 
conclusions in $6] 

We wish to emphasize from the outset that we are in- 
terested in how mean flows, namely differential rotation 
and meridional circulation, respond to turbulent stresses 
in the solar NSSL that are currently unknown. Thus, in 
expressing the relevant dynamical balances we will write 



terms involving the inertia of the mean flows (includ- 
ing the Coriolis force) explicitly. Meanwhile, we will in- 
corporate turbulent stresses by means of two relatively 
generic terms, J- and Q, introduced in sections 12.11 and 
13.11 (see also Appendix A). The reader who values speci- 
ficity may regard T and Q as an embodiment of the con- 
vective Reynolds stress, with explicit expressions given 
in Appendix A. This is likely to be the dominant con- 
tribution in the solar NSSL. However, by leaving these 
terms unspecified we wish to highlight the generality of 
the physical processes we consider. Thus, the Maxwell 
stress and the large-scale Lorentz force in stars and the 
artificial viscous diffusion in numerical simulations can 
establish mean flows in much the same way. 

2. GYROSCOPIC PUMPING 

2.1. Fundamental Concepts 

We begin our investigation of the dynamical balances 
in the NSSL with the standard (compressible) equations 
of magnetohydrodynamics (see Appendix A) , alternating 
throughout between spherical polar coordinates (r, 9, (f>) 
and cylindrical coordinates (A, cj>, z). Thus, A = rsin9 
is the cylindrical radius and z — r cos 9 is the coordinate 
parallel to the rotation axis (each with corresponding 

unit vectors, e.g. A). 

An equation for the conservation of angular momen- 
tum can be readily obtained by multiplying the zonal 
(4>) component of the momentum equation by the mo- 
ment arm, A, and then averaging over longitude and time 
(denoted by angular brackets <>). This yields 



. . dC , . . 2 5S1 . 



■VC + F 



(1) 



where p is the mass density and v is the bulk velocity, 
including meridional and zonal components; v = v m + 
v^cj). The term T incorporates all non-axisymmetric, 
magnetic, and viscous effects, as described further below 
and in Appendix A. Thus, we refer to it as the net axial 
torque. However, if magnetic fields and viscous diffusion 
are neglected, T has only one component, and that is 
the Reynolds stress. We use an inertial (non-rotating) 
coordinate system so the angular velocity is given by f2 = 
(v<t>) A -1 and the specific angular momentum is C = A 2 f2. 
The meridional velocity can be expressed in spherical and 
cylindrical coordinates as v m = v r f + vq6 = v\X + v z z. 

If we assume a statistically stationary state equation 
(fT]) becomes fe.g. lMiesch fc Toomrell2009D : 



(pv m ) -VC = T 



(2) 



Furthermore, if we assume that there is negligible mass 
flux through the solar surface, then a steady solution is 
only possible if 



TdV = 



(3) 



where V is the volume of the entire solar interior r < R. 

We emphasize that equations (Q} and ([2]) are valid for 
any arbitrary value of the Rossby number. They follow 
directly from the MHD equations, as demonstrated in 
Appendix A. The left-hand-side includes all terms involv- 
ing the inertia of the mean flow and the right-hand-side. 
F, is dominated by the Reynolds stress. 
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More specifically, the left-hand-side of equation ^ 
repesents the advection of angular momentum by the 
mean meridional circulation and the right-hand-side in- 
corporates all other forces, most notably the convective 
Reynolds stress and the Lorentz force. Molecular viscos- 
ity can also contribute to the net torque J- but this is ex- 
pected to be negligible in stars. Any imbalance between 
the forces that contribute to T will induce a meridional 
flow (pv m ) across C isosurfaces. Although ft contours are 
nearly radial in the bulk of the solar convection zone, C 
contours are nearly cylindrical, as demonstrated in Fig- 
ure 1, increasing away from the rotation axis such that 



(4) 



and dC/dX > 0. Thus, equation ([2]) implies that a ret- 
rograde torque T < will induce a flow toward the ro- 
tation axis while a prograde torque T > will induce a 
flow away from the rotation axis. This is the concept of 
gyroscopic pumping discussed by Mclntyre (1998, 2007; 
see also Haynes et al. 1991, Garaud & Acevedo Arreguin 
2009, and Garaud & Bodenheimer 2010). Thus, we refer 
to equation (j2J as the gyroscopic pumping equation. 

The process by which equation © is established will 
be addressed in Sj3j Here we focus on its implications. To 
proceed, we assume that the mass flux is divergenceless, 
V- (pv m ) = 0. This can be justified with the anelastic 
approximation but that is not necessary; a non-divergent 
mean mass flux follows directly from the compressible 
mass continuity equation under the assumption of a sta- 
tistically stationary state. We may then define a stream- 
function such that 



<9* 

(PV\) = -TT and (P V z 

oz 
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XdX 



(5) 



If the C contours are cylindrical, as expressed in equa- 
tion (QJ, then follows directly from equations © and 
©: 



*(A,z) 



dC 
dX 



F{X lZ ')dz' 



(6) 



where z b = (i? 2 — A 2 ) 1 / 2 . In obtaining equation © we 
have assumed that there is no mass flux through the 
photosphere so if? = at r = R. Thus, *f? is obtained 
by integrating along a cylindrical surface, beginning at 
the photosphere in the northern hemisphere (z = Zf,) and 
proceeding in the negative z direction, through the equa- 
torial plane, and continuing to the photospheric bound- 
ary in the southern hemisphere (z — —Zb)- 

An appreciation for how gyroscopic pumping operates 
is best obtained by considering the simplest case in which 
the differential rotation is weak. This limit is satisfied 
if R° R « 1, where R° R = Afi/(2fi) is the Rossby 
number associated with the differential rotation, and Afl 
is some measure of the variation of the rotation rate 
with latitude and radius (a rigorous derivation of the 
limit yields Af2 = A|Vf2|). This criterion is approxi- 
mately satistied in the Sun; helioseismic inversions indi- 
cate R® R ~ 0.16. In the low Rossby number limit, the 
uniform rotation component £Iq dominates the angular 
momentum gradient and dC/dX — 2A51o. Substituting 
this into (|6j) then yields *f? for a given T . 



This simple example emphasizes an important point: 
Here, the net zonal force (encompassed in the net ax- 
ial torque J-) determines the meridional flow, not neces- 
sarily the differential rotation. For a given J 7 , equation 
(|2|) provides a direct link between the meridional flow 
and the turbulent angular momentum transport, valid 
to lowest order in the Rossby number (which is assumed 
to be small). The differential rotation is determined by 
other mechanisms, as discussed in £12.21 ££51 an d $5] An 
idealized demonstration of how this works is given in Ap- 
pendix B. 

Note also that a cylindrical torque T = !F(X) is ruled 
out by equation (|6]). This would imply a cylindrical mass 
flux (pv\) that is independent of z, which is in turn ruled 
out by mass conservation and our requirement that there 
be no flow through the surface r — R. More explicitly, we 
can say that if the C profile is cylindrical [eq. (j4}] , then a 

steady meridional flow is only possible if J Zb Zb J~dz — 0. 
In other words, the meridional flow responds mainly to 
the axial variation of the torque, d? ' / 1 dz. The amplitude 
of the resulting v TO is proportional to T and inversely 
proportional to fio! For a given rotation rate, the stronger 
the zonal force, the stronger the meridional flow that is 
induced. 

For finite values of R^ R equation ([2]) must still hold 
in a steady state but the contribution of the differential 
rotation to V£ cannot be neglected. The meridional cir- 
culation will redistribute angular momentum so C will 
depend on "J and the gyroscopic pumping equation is 
nonlinear even if J- is fixed. Linear and nonlinear feed- 
backs of il and $ on J complicate the problem further. 
A unique solution requires consideration of the merid- 
ional momentum and energy equations, as well as mean 
flow profiles that adjust in order to minimimize the net 
axial torque J 7 . We address these issues in £12.21 and £j3] 
below. 

2.2. Clarifications and Generalizations 

As emphasized in the last paragraph of £12.11 the gyro- 
scopic pumping equation ([2]) is robust. It holds for any 
arbitrary value of the Rossby number provided there ex- 
ist well-defined, persistent mean flows. Furthermore, V£ 
is undeniably directed away from the rotation axis in the 
NSSL, as revealed by global helioseismic inversions (Fig. 
Q})). Thus, the sense and amplitude of the meridional 
flow is linked to the net axial torque J- '. However, as 
also noted in £12. 1[ this link is in general nonlinear and 
depends on factors other than the angular momentum 
transport. In this section we discuss some of the subti- 
tles involved. 

We begin with equation (|6]), which rests on the approx- 
imation of a cylindrical angular momentum profile, as ex- 
pressed in equation (j4|). Although this is an instructive 
example with relevance to the Sun (cf. Fig. [TJ)), the con- 
cept of gyroscopic pumping is more general and applies 
straightforwardly to other scenarios as well. In general, 
the meridional flow through C isosurfaces is linked to the 
net axial torque T while the flow along C isosurfaces fol- 
lows from mass conservation. The integral in equation 
([6]) for "J would then proceed along C isosurfaces. 

Another issue mentioned in £12.11 concerns the depen- 
dence of the net axial torque J- on the mean flows and 
the possibility that the mean flows adjust to minimize T . 
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In some circumstances, this dependence may be regarded 
in terms of a differential operator that operates on the 
mean rotation profile: J- = J-{tt}. An example is the 
case of turbulent diffusion, in which the angular momen- 
tum flux is proportional to as discussed in Sj5] In 
this case, the gyroscopic pumping equation @, admits 
a homogeneous solution as well as a particular solution. 
In other words, we may write Q = £l) t + tt p where 



F{Q h } = and T{n p } 



(pv m ) -V£ 



(7) 



This decomposition is only rigorously valid if the oper- 
ator is linear and if the dependence of the angular mo- 
mentum advection (pv m ) C on f2 is linearized in some 
way (e.g. for R® R << 1). Even so, it serves to illustrate 
intuitively how the nature of J- has implications for the 
differential rotation as well as for the meridional circula- 
tion. Even in the absence of meridional flow, the solution 
to the homogeneous equation T = may in general ex- 
hibit a differential rotation that depends on the nature 
of the operator and the boundary conditions. Examples 
are given in ^5.21 

Yet, the basic premise of gyroscopic pumping is still 
valid, namely that a net axial torque J- ^= can only be 
sustained in a steady state if angular momentum is con- 
tinually replenished by advection from the surrounding 
fluid. The meridional flow, differential rotation, and tur- 
bulent stresses will adjust until this is achieved. Merid- 
ional and zonal forces both play a role in this nonlin- 
ear dynamical adjustment as described in £13.21 and con- 
tribute to the mean flow profiles that are ultimately re- 
alized. Still, the nearly cylindrical C profile revealed 
by helioseismology (Fig. QJ) provides a robust link be- 
tween the meridional flow and the net zonal forcing, ex- 
pressed in equation Furthermore, the prominent 
axial £1 gradient dil/dz revealed by helioseismic inver- 
sions provides a robust link between the rotational shear 
and the meridional forcing, as addressed in <JJ] Such 
non-intuitive links between meridional/zonal flows and 
zonal/meridional forcing are mediated by the Coriolis 
force so they are most prominent in rapidly rotating sys- 
tems. However, their robustness applies even in the solar 
NSSL where the Rossby number based on convection is 
large (Q. 

Another important point about gyroscopic pumping 
is that it is inherently a non-local process. A local- 
ized torque T will in general induce a glob al circulation, 
exten ding far beyond the forcing region (jHavnes et al.l 
I1991D . This has particular significance with regard to 
the problem of tachocline confinement, whereby gyro- 
scopic pumping in the convection zone induces a cir- 
culation that burrows downward into the radiative in- 
terior w ith time unless othe r physical processes sup 
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Garaud fc BodenheimerlfeOI Of ). Likewise, zonal forces in 



the NSSL have potential implications for mean merid- 
ional and zonal flows throughout the convection zone. 
For a demonstration of how a local zonal force in the 
NSSL can induce a global meridional flow, see the ana- 
lytic example presented in Appendix B. 

Yet, the coupling between the bulk of the convection 
zone (CZ) and the NSSL will also work the other way. 
Namely, there must be a net torque in the convection 



zone Tcz that induces a meridional flow by gyroscopic 
pumping that will extend into the NSSL. Furthermore, 
given the relatively large mass and energy content of 
the CZ relative to the NSSL, we may expect this merid- 
ional flow to overwhelm that which is driven by net axial 
torques within the NSSL, Tnssl- In order to account 
for the solar differential rotation, the sense of the deep- 
seated angular momentum transport must be such that 
Tcz is positive at low latitudes and negative at high lat- 
itudes. This will induce a counter-clockwise circulation 
in the northern hemisphere via eq. @ that will pump 
mass flux into the NSSL, maintaining a poleward flow 
even in the absence of any turbulent stresses within the 
NSSL itself. However, if Tnssl were indeed zero, then 
this CZ circulation would redistribute angular momen- 
tum until the C contours aligned with the streamlines 
of the meridional flow such that (pv m ) -V£ = 0. This 
is clearly not the case in the Sun, where the meridional 
flow is poleward and the C contours are nearly cylin- 
drical (Fig. [TJ)). The meridional flow clearly crosses C 
isosurfaces, so Tnssl must be nonzero. 

In short, gyroscopic pumping by net axial torques in 
the deep convection zone may contribute to the pole- 
ward flow in the NSSL but angular momentum trans- 
port within the NSSL itself must also play a ro le. This 
is con sistent with the mean- field simulations of iRempel 
(2005) who considered the maintenance of mean flows 
in the solar interior based on idealized parameteriza- 
tions for the turbulent momentum and energy transport. 
He found poleward meridional flow near the solar sur- 
face even without a NSSL. When he included an inward 
(cf. i)5.4[) angular momentum transport in a thin layer 
near the surface, he found both an NSSL-like shear layer 
where dtt jdr < as well as an enhancement of the pole- 
ward flow. 

2.3. Implications 

In light of the discussion in $'2.1\ and ^2.21 the implica- 
tions of equation ^ for the solar NSSL are clear. The 
same physical mechanism responsible for the deceleration 
of the rotation rate in the solar surface layers inferred 
from helioseismology (dfl/dr < 0) is also responsible for 
the poleward meridional flow inferred from helioseismic 
and Doppler measurements ((vg) < in the northern 
hemisphere). The underlying cause of both phenomena 
is a retrograde net axial torque (zonal force) T which 
is most likely due to a divergence in the angular mo- 
mentum transport by the convective Reynolds stress. Its 
localization near the solar surface (resulting in a large 
dT/dz that effectively generates a meridional flow) is in 
turn a likely consequence of the rapidly changing length 
and time scales of convection, L c and t c . The influ- 
ence of rotation on convection is typically quantified by 
the Rossby number based on the convective time scales 
R = (2r2r c ) _1 (as opposed to R^ R above). Estimates 
for giant cells (r ~ 10-20 days, R ~ 0.1-0.2) and gran- 
ulation (r ~ 8 min, R Q ~ 400) suggest that R should 
cross unity somewhere in the vicinity of the lower NSSL, 
likely signifying a qualitative change in convective trans- 
port. 

Thus, in this section we have established that a net 
retrograde zonal force J- < in the NSSL will induce a 
poleward meridional flow. This is a very general result, 
independent of the nature of J-, which we address in 
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Before proceeding to this, however, we discuss another 
general and important maxim; a net axial torque T can- 
not in itself account for the existence of the NSSL. Again, 
this conclusion is completely independent of the nature 
of T and follows essentially from the Taylor-Proudman 
theorem. Coriolis-induced meridional flows will tend to 
eliminate all axial rotational shear dQ/dz regardless of 
the amplitude and structure of the zonal forcing. Thus, 
the axial shear dfl/dz in the NSSL must be determined 
not by convective angular momentum transport T but 
rather by turbulent and baroclinic stresses in the merid- 
ional plane, which we now address (jjBJ). 

3. MERIDIONAL FORCE BALANCE 

3.1. The Zonal Vorticity Equation 

In i 32.ll we discussed how a net axial torque T can 
induce a meridional circulation through the gyroscopic 
pumping equation (|2|) but we said little about how this 
balance is achieved. We also left open the question of how 
the differential rotation profile is established, particularly 
in the low Rossby number limit (R® R << 1) when T is 
directly linked to the meridional flow through eq. ((6|). 

To address these issues for any arbitrary value of the 
Rossby number, we must also consider the meridional 
components of the momentum equation. Exploiting the 
divergenceless nature of the mean mass flux (pv m ), we 
can combine these two equations into one by considering 
the zonal component of the curl; in other words, the zonal 
vorticity equation, averaged over longitude and time. In 
order to illustrate how the system adjusts in response to 
specified forcing scenarios, we will temporarily retain the 
time derivative and write this equation as follows: 



■d_ 

Of 



dfl 2 

(w*) = A— — + B + Q 
oz 



(8) 



The vorticity is defined as w = Vxv, with uj^ as the 
zonal component. We emphasize that equation ((8} is 
valid for any arbitrary value of the Rossby number, as 
demonstrated in Appendix A. The first term on the right- 
hand side includes the complete inertia of the mean zonal 
flow, including uniform (Coriolis) and differential rota- 
tion components. This is followed by the baroclinic term 



B : 



V(p)xV(P) 

(P) 2 







(9) 



where p is the density as before and P is the pressure. 
The final term on the right-hand-side of equation flSJ), 
g, represents turbulent stresses in the meridional plane. 
As with the net axial torque J 7 , Q includes the Reynolds 
stress, the Lorentz force, and the viscous diffusion. We 
have also included in Q baroclinic contributions associ- 
ated with thermal fluctuations, e.g. p' = p — (p). An 
explicit expression for Q is given in Appendix A, along 
with the derivation of equation ([5]). 

3.2. A Thought Experiment 

We now describe a simple thought experiment in order 
to illustrate several important points about gyroscopic 
pumping and to gain an intuitive feel for how it operates 
within the context of the NSSL. We stress that this is 
an idealized example intended to illustrate fundamental 
physical principles, namely that some meridional forc- 
ing is needed to account for the existence of the NSSL. 



We then proceed to discuss the nature of this meridional 
forcing in §3.3. 

Consider a rotating spherical volume of radius R sub- 
ject to a specified net axial torque T that turns on at 
some instant t = to. In analogy to the NSSL, we will 
assume that T < in a thin layer near the surface, say 
r s < r < R. To be explicit, we can set r s 0.95i? 
(see £|5.2[) . In order to allow the system to eventually 
reach an equilibrium state, we assume that the net (in- 
tegrated) torque in the convection zone is positive, such 
that T satisfies equation |3]). 

For simplicity we will assume that the baroclinic and 
turbulent stresses in the meridional plane vanish, so 
B = Q = 0. This corresponds to a flow that is ax- 
isymmetric, non-magnetic, non- diffusive, and adiabatic. 
Furthermore, we assume that J- = for t < to- Thus, 
before time to, the system can sustain an initial equilib- 
rium state with a cylindrical rotation profile = f2i(A) 
and no meridional flow, v m = (satisfying the Taylor- 
Proudman theorem; see, e.g. Pedlosky 1987). In analogy 
with the Sun, we assume that the equator rotates faster 
than the poles, so dili/dX > 0. It follows also that the 
angular momentum increases outward, V£-A > and 
is cylindrical (V£-z = 0), as expressed in equation Q. 
This initial equilibrium satisfies equations and @, 
with d (cjj,) jdt — 0. 

Now turn on the torque T at t = to. How does the 
system respond? Initially, the torque will slow down 
the rotation rate in the NSSL, tilting the f2 contours 
away from the rotation axis. This will produce an ax- 
ial rotation gradient dfl/dz which will in turn induce a 
counter-clockwise meridional circulation in the northern 
hemisphere (negative (w^)) according to equation @. In 
other words, there will be a poleward flow in the NSSL 
and a net equatorward return flow in the deeper convec- 
tion zone to ensure mass conservation. 

This induced meridional flow will redistribute angular 
momentum, altering the profile. In the NSSL, angular 
momentum will be advected from low latitudes toward 
the poles, accelerating the rotation rate, accompanied by 
a deceleration of the deeper convection zone. This will 
proceed until a new, final equilibrium is reached, again 
with a cylindrical rotation profile Q = flf(X). The new 
profile will be different from the initial profile (f2/ ^ fii), 
with cylindrical isorotation surfaces shifted away from 
the rotation axis (positive A direction). Furthermore, 
the torque J- will sustain a steady, nonzero meridional 
circulation \& f , given by equation ^ , with poleward flow 
in the NSSL and an equatorward return flow in the deep 
CZ. 

In summary, the response of the system to a retro- 
grade zonal torque J- localized in the NSSL is to estab- 
lish a poleward flow that will increase in amplitude until 
the prograde angular momentum advected into the layer 
from the deep CZ balances the retrograde angular mo- 
mentum imparted by the torque. The meridional flow is 
predominantly poleward (v m w VgO) rather than cylin- 

drically inward {v\\) because of the thin radial extent of 
the NSSL (see Appendix B for a demonstration) . Mean- 
while, the rotation profile is altered but still cylindrical, 
n = 0/(A). 

This simple thought experiment illustrates two impor- 
tant points. First, gyroscopic pumping is mediated by the 
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Coriolis force but the two phenomena are not equivalent. 
The meridional components of the Coriolis force can in 
principle vanish (for the special case of flf = constant) 
while the system still sustains a gyroscopically-pumped 
meridional circulation. More generally, for any cylin- 
drical rotation profile fi(A), the meridional components 
of the Coriolis force are opposed by pressure gradients 
(geostrophic balance) so they impart no zonal vorticity 
and thus no meridional momentum. The meridional flow 
is determined not by the Coriolis force (which responds 
to O), but rather by the conservation of angular momen- 
tum (which responds to J-). 

The second point is even more important: The mere 
existence of a retrograde torque J- does not guarantee 
the presence of a near-surface shear layer. As stated 
in the introduction, even if the concept is valid, angular 
momentum conservation by turbulent convection alone 
cannot account for the existence of the NSSL. In the 
absence of baroclinic, turbulent, or magnetic stresses in 
the meridional plane (B = G = 0), the meridional circu- 
lation will wipe out all axial shear dfl/dz, regardless of 
the nature, magnitude or profile of the connective angular 
momentum transport (embodied by J-). 

The robustness of cylindrical rotation profile s is a con- 
seque nce of the Taylor-Proudman theorem (e.g. lPedloskvl 
1987). No matter how strong T is, a commensurate 
meridional circulation will suppress axial shear on a time 
scale comparable to the rotation period (~ 28 days). 
The prominent axial shear dfl / dz exhibited by helioseis- 
mic rotation inversions (Fig. [T]i) then implies that some 
meridional forcing, either baroclinic or "turbulent" , B or 
G, is necessary to account for the existence of the NSSL 
as well as the detailed structure of the inferred ft profile. 

3.3. What Determines the Rotational Shear in the 
NSSL? 

In this section we argue that the existence and location 
of the NSSL can be attributed to a qualitative change in 
the meridional force balance, marked by a transition from 
baroclinic to turbulent stresses; B to Q. Furthermore, 
it is the meridional stresses (G and/or B) rather than 
the angular momentum transport that largly determine 
the slope of the profile (although the negative sign of 
dQ/dr still requires a retrograde net axial torque F < 0). 

Although there are still some subtle uncertainties re- 
garding the underying dynamics, recent global solar con- 
vection simulations and mean-field models have con- 
verged on a consistent paradigm whereby the nearly ra- 
dial angular velocity contours (conical 57 isosurfaces) at 
mid-latitudes in the deep CZ are attributed to baroclinic 
forcing 



on 2 

dz 



B 

A 



9 9(S) 
rXCp 89 



(deep CZ) (10) 



where g is the gravitational acceleration, S is the specific 
entro py, and Cp is the specific he a t at constant pres- 
sure (IKitchatinov fc Riidigerl 119951: lElliott et all 2000; 



Robinson & Chan 2001; Brun fc Toomre 2002; Rempcl 
20051 : iMiesch et al.l 120061 : iBrun et all 1201 lft. This is 



referred to as thermal wind balance (IPedloskvl fl^87t 
IMiesch fc~T oomrc 2009) . The convective Reynolds stress 
is still necessary to account for the amplitude and 
sense of the angular velocity contrast between equator 



and pole, Af2, but without baroclinicity, models gen- 
erally yield cylindrical (Taylor-Proudman) profiles, in 
contrast to the conical profiles inferred from helioseis- 
mic inversions (Fig. [T]i). In early models, the latitudi- 
nal entropy gradient in equation (|10[) is established en- 
tirely by a latitude-dependent convective heat flux at- 
tribut ed to the influence of rotation on convective mo- 
tions (IKitchatinov fc Riidigerl 119951: lElliott et al.l 120001: 
iRobinson fc Chan! 120011 : IBrun fc Toomrd l2002( l How- 
ever, recent models have demonstrated that thermal cou- 
pling to the subadiabatic portion of the tachocline by 
means of a gyroscopically-pumped meridional circula- 
tion may c ontribute to establishing the requisite thermal 
gradi ents (|Rempell 120051: IMiesch et aT] 120061: IBrun et al.l 
[20III). 

In any case, the nature and prominence of the baro- 
clinic term in convection simulations and mean-field 
models hinges on the low effective Rossby number in the 
deep CZ. By effective Rossby number we are referring to 
the amplitude of the meridional components of the con- 
vective Reynolds stress relative to the deflection of the 
differential rotation by the Coriolis force and the baro- 
clinic forcing. However, as noted in §2.3| the Rossby 
number rises steadily with radius in the convection zone, 
becoming much greater than unity in the NSSL. This 
alone suggests that the turbulent stresses, represented 
by G, may overwhelm the baroclinic term B in the solar 
surface layers. In fact, this shift in the meridional force 
balance (from B to G) may mark the base of the NSSL 
even more significantly than a change in the turbulent 
angular momentum transport J- . 

Thus, we propose that in the NSSL thermal wind bal- 
ance, equation (fTU|) is replaced by 



on 2 

dz 



G 
A 



(NSSL). 



(11) 



The implication would then be that turbulent stresses 
in the meridional plane, G, largely determine the dif- 
ferential rotation profile in the NSSL. This hypothesis 
is consistent with global simulations of solar convection 
which, although they do not yet reproduce the NSSL, 
do exhibit a departure from therm al wind balance in 
the surface layers abo ve r ~ 0.95 (lElliott et al.l 120001 : 
IBrun et al.ll2010l 120111 ). This departure is attributed to 
a combination of resolved convective Reynolds stresses 
and modeled sub-grid scale diffusion, both captured in 
our turbulent transport term G- Achieving a NSSL in 
global convection simulations will likely require a reli- 
able representation of turbulent transport by unresolved 
(subgrid-scale) motions, as discussed in $2] and $5] 

Can solar observations and models reveal what the 
meridional force balance is in the NSSL, or in other 
words, which is bigger, B or CJ? In principle the answer 
is yes, but as we discuss in the remainder of this section, 
current results are inconclusive. 

Helioseismic rotational inversions provide a good es- 
timate for the left-hand-side of equation (p~0|) . namely 
C = —\dil 2 /dz (Q. If we had a reliable estimate for 
the right-hand-side, B, from solar observations and if its 
amplitude were << C, then we could rule out equation 
(fT0|) and thereby verify the proposed balance expressed 
in equation (jlip . 

One way to estimate the amplitude of B, at least in 
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principle, is by means of helioseismic structure inver- 
sions. The thermal gradients that give rise to B imply 
aspherical sound speed variations. If these gradients are 
to balance the axial shear in the rotation rate through 
baroclinic torques as expressed in equation (fTUj) . then 
their relative amplitude must be of order 



AS 



2 d (S) 2r\ dn 2 



itCp 86 



*g 



dz 



2 x 10 



-5 



(12) 



where AS is the entropy variation from equator to pole. 
The numerical estimate in equation (|13|) is based on 
global rotational inversions as described in SJH below. 

Such subtle thermal variations are beyond the de- 
tection limit of current heliosei smic structu r e in- 
versions, even m t he N SS L (IGoueh et alJ [1996; 
IChristensen-Dalsgaardl I2002D . IBrun et all (|2010D have 
reported aspherical sound speed variations much larger 
than this, of order 10 _5 -10 -4 , but the amplitude of the 
variations peaks near the base of the convection zone at 
high latitudes, where the reliability of the inversions is 
questionable. If these inversions are indeed valid and 
if they do indeed trace thermal gradients as opposed 
to magnetic effects, then thermal wind balance, equa- 
tion (ITU)) would not be satisfied and turbulent stresses Q 
would have to contribute to the meridional force balance 
throughout the deep convection zone as well as the NSSL. 
More work is needed to set stricter limits on aspherical 
sound speed variations from helioseismic structure inver- 
sions and whether they are consistent with helioseismic 
rotational inversions within the context of thermal wind 
balance. For now we regard the issue as not yet settled. 

An independent estimate for B near the solar surface 
can (again, in principle) be obtained from photospheric 
irradian.ee measurements. The final equality in equa- 
tion (JTUJ, expressing B in terms of the latitudinal en- 
tropy gradient, applies for an ideal gas equation of state 
and a nearly hydrostat ic, adiabatic stratification (e.g. 
iMiesch fc Toomrel 12009ft . An alternative expression can 
be obtained by assuming that the photosphere is an iso- 
baric surface. In this case the baroclinic term can be 
expressed in terms of the photospheric temperature vari- 
ation from equator to pole, AT. If this is to balance the 
Coriolis term C, then we must have relative temperature 
variations of the same order as in equation (|12[) : 



AT 2d_ 2RX dtt 2 

T ~ 7T 06 n ^ ' ~ irg dz 



2 x 10" 5 



(13) 



We have again assumed an ideal gas equation of state and 
hydrostatic balance for simplicity. If we further assume 
that the solar irradiance arises from blackbody emission, 
then we should expect to see a corresponding latitudinal 
variation of the irradiance. In other words, if the solar 
differential rotation in the NSSL were indeed in thermal 
wind balance [eq. (|10j) ]. then we would expect the poles 
to be brighter than the equator with a relative irradiance 
enhancement of AI/I ~ 8 x 10~ 4 . This corresponds to 
a temperature difference of about 0.1K. 

Irradiance variations of a t least this order have indeed 
been detected (reviewed bv lRast et al.ll2008ft. In partic- 
ular, recent measurments bv iRast et al.1 (|2008ft indicate 
that the poles are indeed brighter than the equator, with 
AI/I ~ 1.5 x 10 -3 . However, whether or not these can 
be interpreted as temperature variations is questionable. 



There is a strong possibility that the irradiance variation 
arises instead from unresolved magnetic flux elements. 
Thus, we view current estimates of B based on photo- 
spheric irradiance measurements as inconclusive as well. 

The essential question is whether B can keep up with 
the sharp increase in dVl/dz throughout the NSSL. The 
answer depends on the details of convective heat trans- 
port in the solar surface layers, which are not well under- 
stood. If convection acts as a turbulent thermal diffusion, 
then the thin shell geometry implies that the latitudi- 
nal temperature gradient (AT) at the base of the NSSL, 
r = r s , might be efficiently transmitted to the surface. 
The rapid decrease of the background temperature with 
radius might then lead to an enhancement of the relative 
latitudinal temperature gradient AT/T, which in turn 
would imply an increase in B [cf. eq. (|13[) ]. However, 
the horizontal heat transport and poleward meridional 
flow would tend to suppress latitudinal thermal gradi- 
ents (AT) and thus diminish B. 

To approach this issue more quantitatively, we esti- 
mate from helioseisimic inversions (SJU) that C is a fac- 
tor of 5-10 larger at the surface than at the base of the 
NSSL, r = r s ~ 0.946i? (defined as where the mean ra- 
dial f2 gradient changes sign; see §5.2| . Over that same 
radial range the background temperature decreases by 
a factor of 30-40. Thus, taking into account also the 
variation of g and r, we estimate that horizontal trans- 
port would have to be at least three times more efficient 
than vertical transport in order for the magnitude of B 
to increase less rapidly than C, assuming thermal wind 
balance at r — r s and that the photosphere is an iso- 
baric surface, as in equation (|13[) . We may expect the 
ratio of horizontal to vertical turbulent diffusion to scale 
as Kh/n v ~ UhLh/(U v L v ) where U and L are typical 
velocity and length scales and h and v denote horizonal 
and vertical directions. Mass conservation, V-(pv) = 
implies Uh/U v ~ L^/ 'L v , so Kh/ 'k v ~ (Lh/L v ) 2 . It is rea- 
sonable to expect that this ratio may indeed exceed three 
for solar surface convection. For example, if we put in 
numbers for granulation, ~ 1 Mm and L v ~ H p ~ 300 
km, where H p is the density scale height, then Kh/ K v ~ 
10. This helps to substantiate the conjecture that baro- 
clinic forces may be too weak to balance the Coriolis force 
associated with the differential rotation in the NSSL. 

In summary, we suggest that the characteristic change 
in the slope of the profile that defines the NSSL is 
largely determined by a shift in the meridional force bal- 
ance from baroclinic to turbulent stresses; that is, from 
equation (fT0|) to equation (fTTj) . Although this hypothe- 
sis is in principle testable, we believe that current data is 
insufficient to either confirm or deny it. In the remainder 
of this paper we assume that this transition does indeed 
occur and we explore its implications with regard to the 
dynamics of the NSSL. 

Note also that the statements made in this section 
apply mainly to the persistent, background component 
of the mean flow that is present throughout the solar 
cycle. The baroclinic term B may well establish time- 
varying meridional and zonal flows such as the low- 
latitude br anch o f the so lar to rsional o s cillati on as dis- 
cussed by iSpruitl (|2003t) and iRempell (|2007ft . In the 
Spruit-Rempel model, enhanced cooling in magnetically 
active latitudinal bands induces converging meridional 
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flows that in turn accelerate zonal flows by means of the 
Coriolis force. If the cooling is sustained for at least sev- 
eral months (long enough for the Coriolis force to fully 
respond), then the dynamical balances in equation ([2]) 
and ([TUl) can be established in a quasi-static sense. As 
for the persistent background flow, the coupled inertia of 
the fluctuating flow components provides a tight link be- 
tween the meridional circulation and the rotational shear. 

4. HELIOSEISMIC ESTIMATES OF TURBULENT 
TRANSPORT 

In <21 and <J3] we argued that the dynamical balances 
that are likely to prevail in the NSSL are given by equa- 
tions ^ and (fTTj) . The left-hand side of each of these 
equations involves quantities that are accessible to he- 
lioseismic inversions, namely Q and (v m ). Thus, if we 
assume these balances hold, then we can obtain observa- 
tional estimates for the turbulent stresses T and Q. 

That is the focus of this section. We will obtain es- 
timates for the turbulent stresses F and Q based on lo- 
cal and global helioseismic inversions and we will discuss 
their implications for the nature of turbulent transport in 
the NSSL. Then, in $5] we will investigate whether these 
results are consistent with simple theoretical models. We 
begin in £14.11 by discussing the data we use and we then 
proceed to calculate T and Q in 



4.1. Helioseismic Data and Inversions 

We use both local and global helioseismic procedures to 
estimate the solar rotation rate f2 and the meridional flow 
component (v m ). Both techniques exploit the fact that, 
through the Doppler effect, a flow induces a frequency 
splitting between acoustic waves propagating in opposite 
directions. By measuring this frequency splitting for a 
variety of waves that reside within different regions of the 
solar interior, through an inversion process, maps of the 
flow as a function of latitude and depth can be produced. 
Global helioseismology employs modes of long horizontal 
wavelength that comprise the sun's global acoustic res- 
onances. The frequency splittings of such modes can be 
used to measure the axisymmetric component of the ro- 
tation rate over a broad range of latitudes (within 70° of 
the equator) and to depths spanning the entire convec- 
tion zone. Local helioseismic techniques, which measure 
short- wavelength waves, do not permit sampling of the 
flows to such high latitudes or to such large depths (only 
within 50° of the equator and only within the upper 15 
Mm or 2% by fractional radius); however, local helioseis- 
mology is capable of measuring both the rotation rate 
and the meridional component of the flow, unlike global 
helioseismology. 

We utilize a local helioseismic procedure called ring 
analysis which assesses the speed and direction of sub- 
surface horizontal flows by measuring the a dvection of 
ambient acoustic waves by those flows (e.g. EH [1988; 
IHaber et~a l. 2002) . Measurement of the frequency split- 
ting provides a direct measure of the fluid's flow velocity 
in those layers where the waves have significant ampli- 
tude. The frequency shift induced by the_flow on any 
given wave component is given by Aw =_k-U, where k is 
the wave's horizontal wavenumber and U is the integral 
over depth of the horizontal flow velocity weighted by 
a kernel which is approximately the kinetic energy den- 
sity of the acoustic wave. The frequency splittings for a 
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Fig. 2. — Rotation profile Q obtained from local and global inver- 
sions, plotted verses (a) latitude (at the indicated radius) and (6) 
radius (at the indicated latitudes). Inversions for the year 1996 are 
shown as a black line and other years 1997-2003 are shown as blue 
symbols. The green curve represents the global rotation inversions 
shown in Figure [l] and the dashed red curve represents an analytic 
fit t o the local 1996 inversions of the form expressed in equation 
JTB1 below. 




-20 
-30 



r/R = 0.990 



30 






20 






10 

















10 






20 


r/R - 0.979 




30 







-40 -20 20 
Latitude (deg) 



40 



-40 -20 20 
Latitude (deg) 



40 



Fig. 3. — Longitudinally-averaged (co-) latitudinal velocity (vg) 
obtained from local helioseismic inversions (ring analysis), plotted 
versus latitude at (a) r = 0.990.R and (b) r = 0.979H. Black lines 
and symbols denote data from 1996 and other colors denote other 
years, as indicated. We focus on 1996 but we show other years 
merely to highlight the magnitude of the variations; the relation- 
ship between meridional flow variatio ns and solar activity has been 
inves t igated in detail elsewhere (e.g. Hathaway 1996; Hab er et al.l 
2002; Zhao & Kosovichcv 2004; Gonzalez-Hernandez et al. 2006; 
Ulrich 2010; Hathawav & Rightmirc 2010, 2011). The red dashed 
line s represent a fit to the 1996 inversions as expressed in equation 
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large number of waves of different wavenumber and ra- 
dial mode order form a system of integral equations that 
can be invert ed to obtain the horizontal flow as a func - 
tion of depth (jThompson et al.lll996t [Haber et al.ll2004D . 
By repeating the analysis over many different locations 
on the solar surface and over many different days, a map 
of the horizontal flow as a function of longitude, latitude, 
depth and time can be generated. 

Our local helioseismic inferences were made using 
Dopplergram data from the Michelson Doppler Im- 
ager (MDI) aboard the Solar and Heliospheric Observer 
(SoHO). Flows have been obtained using the Dynamics 
Campaign data, of which 2-3 months of data per year is 
available. We have chosen to use data from the rising 
phase of the solar cycle from 1996 through 2003. The 
short-wavelength waves used by local helioseismic proce- 
dures are trapped near the solar surface, and the ring- 
analysis measurements used here are therefore capable of 
sampling down to a depth of 15 Mm below the photo- 
sphere, or the upper half of the NSSL. Each local analysis 
samples a region of Sun that is roughly 180 Mm (15° in 
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heliographic angle) in diameter. On any given day, a 
mosaic of 189 overlapping analysis regions are indepen- 
dently analyzed, resulting in a flow map with a horizontal 
spacing of 7.5° between measurement regions, each with 
a horizontal resolution of 15°. 

Estimates for the mean rotation rate SI and the latitu- 
dinal component of the meridional flow (vg), have been 
obtained by averaging the resultant maps over all longi- 
tudes and over the 2-3 months that the Dynamics Cam- 
paign data is available each year. The mean rotation rate 
SI is illustrated in Figure EJ The black curve corresponds 
to measurements made in 1996, while the blue diamonds 
show the weak variability in the flows during the years 
1997 to 2003. The latitudinal component of the flow, 
(vg), is shown in Figure [3] at two different depths and 
with a differently-colored curve for each year. 

The global helioseismic assessments were made using 
four non-overlappi ng 60 day intervals of GONG dat a 
from the year 1996 (H owe et al.ll2000HSchou et al.ll2002f ). 
The rotation rate is inferred from the frequency splitting 
between modes of equal but opposite azimuthal order 
±m. Such modes are identical in all respects other than 
the fact that they propagate in opposite directions in 
longitude around the sun. Solar rotation Doppler shifts 
the mode frequencies in opposite directions and the fre- 
quency splittings of many different modes can be inverted 
to obtain the rotation rate as a function of latitude and 
depth. The result of a regularized Least Squares (RLS) 
inversion is shown in Figures [T] and [2] 

4.2. Estimating J- and Q 

If we assume the dynamical balances in equations ([2]) 
and (jlljl hold, then we can turn them around and cal- 
culate J- and Q based on helioseismic inversions. In par- 
ticular, if we assume that to lowest order p w (p), then 
equation @ yields an estimate of the specific torque 



T = 



T 



dC 
Or 



(vg) dC 
r 39 



(14) 



Thus, the right-hand-side only depends on C = A (v^) 
A 2 S1 and (v m ). Similarly, equation ((TTj) gives 



Q = —X 



dfl 2 

dz 



(15) 



which only depends on SI. 

Estimates for F and Q based on equations (TT4")) and 
(|15|) are shown in Figure 2] In the remainder of this sec- 
tion we provide further details on how these results were 
obtained. For both F and Q we focus on data from 1996. 
Since this corresponds to solar minimum, it minimizes 
the influence of the cyclic component of the solar mag- 
netic field, focusing instead on the persistent turbulent 
stresses that maintain the NSSL throughout the solar 
cycle (E3J|. 

Before proceeding to the more involved calculation of 
T, we first consider the estimate for Q shown in Figure©. 
Since this only depends on fi, we use the global inversions 
shown in Figure [T] to compute it. These extend to higher 
latitudes and lower radii than the local inversions and 
the longitudinal coverage provided by global modes may 
more accurately reflect the axisymmetric SI gradients in 
the NSSL. In order to compute the latitudinal S7 gradient, 
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FlG. 4. — Estimates for the turbu lent stres ses (a) T = Tp^ 1 
and (b) Q, computed from equations ( 1 141) and (1 1 5 I t using local and 
global helioseismic inversions from 1996, as described in the text. 
In (a), the color table ranges from -3.95xl0 8 cm 2 s — 1 (white) to 
zero (black). Blue contours denote negative values and the white 
dotted line indicates the zero contour. In (b) the color scale runs 
from zero (black) to 1.28xl0 -11 s -2 , with blue contours denoting 
positive values. Note that the vertical and horizontal axes are 
different in frames (a) and (6) due to the limited extent of the 
local inversions. The region spanned in (a) is indicated in (b) by a 
white dashed line. A second dashed line indicates the base of the 
NSSL, r s = 0.946-R, defined as where the spherically-averaged Q 
gradient changes sign. 



we first fit the S7 profile to a functional form given by 

n(r, 6) = Sl e (r) + Q 2 (r) cos 2 9 + 4 (r) cos 4 9 (16) 

and then compute the 9 gradient analytically. We then 
compute the radial SI gradient by means of a second-order 
finite difference scheme (on a non-uniform grid) applied 
to the fitting coefficients f2 e , £^2, and SI4. 

The estimation of T from equation (|15[) involves corre- 
lations between the mean meridional and zonal flow. For 
these we use the local inversions described in H4.ll In par- 
ticular, we fit the local SI inversion to a functional form 
as in equation (|16p and the local (vg) inversion (both 
from 1996) to a functional form given by 



(vg) — sin 9 cos 9 (ci(r) + c 3 (r) cos 2 9) 



(17) 



The fits to the local SI and (vg) inversions arc plotted in 
Figures [3] and [2] as red dashed lines. The global SI fit lies 
almost on top of these other curves and is thus omitted 
for clarity. 

The mean radial velocity (v r ) is obtained by assum- 
ing mass conservation; V- (pv m ) — 0. We also assume 
that p ~ (p) ~ p s (r) where p s (r) is the spherically- 
symmetric mean density given by solar structure Model 
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Fig. 5. — Mean radial velocity (v r ) inferred from mass conserva- 
tion, based on local helioseismic inversions for (vg). Solid, dashed, 
and dot-dashed lines correspond to r = 0.977-R, 0.990-R, and 0.996-R 
respectively. Dotted lines indicate zero velocity and the range in 
latitude spanned by the local inversions (± 52.5°). Note that the 
curves extend beyond this range bec ause they are based on the 
analytic expression given in equation J X 9 1) - 



S of IChristense n-Dalsg aard et al.l (|1996l ). This gives 
d , o , ,s rp s d 



di 



(r 2 p s (v r )) 



de 



(sin0 (w e )) 



(18) 



The 9 derivative on the right- hand- side of equation (|18|) 
is computed analytically from the fit in equation (|17[) . 
This leads to the following form for (v r ) 



(v r ) = w e (r) + W2{r) cos 2 9 + Wi(r) cos 4 9 



(19) 



We then proceed by expressing the radial derivative 
on the right-hand-side of equation (fT5]l by means of a 
second-order finite difference scheme and we solve the 
resulting matrix equations for r 2 p s u>i (i — e, 2, 4), with 
boundary conditions such that (v r ) = at r = R. The 
Wi are then smoothed using a boxcar window, taking into 
account the non-uniform grid. The resulting (v r ) profile 
is shown in Figure [5] for several radii. 

Note that the amplitude of the inferred radial veloc- 
ity is extremely small in the NSSL, less than 1 m s _1 
and its contribution to T (and T) is negligible. It does 
contribute a small positive maximum at the equator 
(~ 0.07 x 10 8 cm 2 s" 1 ) as seen in Figure fflx, but the 
prominent maximum at mid-latitudes is due almost en- 
tirely to the final term in equation (I14[) , proportional to 
(vg) dC/d9. Although there is some uncertainty in the 
estimation of (v r ) f ij4.4[) . we believe this conclusion is 
robust. In order for the radial component to make a sig- 
nificant contribution to J 7 , the radial velocity would have 
to be at least an order of magnitude larger than our esti- 
mate. This would imply either a latitudinal velocity that 
is an order of magnitude larger, of order 200 m s , or a 
latitudinal gradient of (vq) that is an order of magnitude 
larger. If the amplitude of the meridional flow is to be 
no larger than 20 ms" 1 , the requisite gradient would re- 
quire that the peak flow be achieved within 3 degrees of 
the equator. Both scenarios can be ruled out by surface 
observations and helioseismic inversions. 

Both T and Q peak at mid-latitudes (Fig. [4} because 
this is where dSl/dz and (vg) dC/89 peak. In fact, these 
two terms (dVl/dz and (vg) dC/89) arc identically zero at 
the equator because of the symmetry implied by equa- 
tions (HHJ) and (fl"7| . As emphasized in i)5.4[ this does 
not mean that the local angular and meridional momen- 



tum flux must vanish at the equator, nor must the mean 
shear. Rather, it implies that the turbulent transport 
and the mean flows adjust themselves in such a way as 
to minimize the net zonal flux divergence (axial torque) 
J- and the net meridional "force curl" Q . This is in con- 
trast to mid-latitudes, where the net zonal and merid- 
ional momentum transport represented by J- and Q must 
be nonzero in order to balance the inertial forces associ- 
ated with the mean differential rotation and meridional 
circulation. This delicate, nonlinear, nonlocal interplay 
between turbulent transport and mean flows ultimately 
determines the structure of the NSSL. 

Note that the sign of F (and T) is negative at mid- 
latitudes, signifying that turbulent stresses must exert a 
retrograde net axial torque, removing the angular mo- 
mentum that is imparted by the meridional flow. The 
sign of Q, meanwhile, is positive, tending to induce a 
clockwise circulation in the northern hemisphere in or- 
der to offset the counter-clockwise circulation that would 
tend to establish a cylindrical (Taylor-Proudman) rota- 
tion profile (gSH). 

4.3. Anisotropy and Energetics 

The quantities plotted in Figure SI T = J ' p~ x and £/, 
represent turbulent (non-axisymmetric) zonal and merid- 
ional momentum transport but a direct comparison be- 
tween them is somewhat elusive since they apply in dif- 
ferent contexts; J- appears in the angular momentum 
equation while Q appears in the zonal vorticity equation; 
their units are different. In this section we put T and Q 
on an equal footing by converting them to accelerations, 
with units of cm s~ 2 . In other words, we can write the 
momentum equation as 



6>< 



dt 



= -«v) -V) (v) 



(p> 



(20) 



Note that equation (|20[) is exact in the sense that it 
is faithful to the compressible MHD equations, with all 
non-axisymmetric influences (as well as some axisymmet- 
ric influences such as viscous diffusion and the mean-field 
Lorentz force) encapsulated in the "turbulent" accelera- 
tion term A. 

Now consider again our dynamical balance equations 
([2]) and (fTTj) . The left- hand-side of both of these equa- 
tions arises from the first term on the right-hand-side of 
equation (|20l) . namely the inertia associated with mean 
flows (we use an inertial reference frame so the Coriolis 
force is included in this term as well; see Appx. A). The 
terms on the right-hand-side of equations ([2]) and (|TT1) . 
T and G, arise from A so we can compute A from the 
helioseismic results presented in Figure SJ Note that the 
second term on the right-hand-side of equation (f2"U|) gives 
rise to the baroclinic term B in equations (JSJ and (ITU| . 
which we propose is negligible in the NSSL fi j3.3[) . 

The zonal component of A is directly related to T: 



.4, 



T_ 



(21) 



In order to obtain a relationship between Q and the 
meridional components of A, A m , we note that any ar- 
bitrary, axisymmetric vector field with only meridional 
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TABLE 1 

Energetics and Time Scales' 1 
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0.985 



0.980 



20 

latitude 



Fig. 6. — Acceleration terms (top) A r , (middle) Ag, and (bottom) 
—Afi, for the region spanned by the local inversions (r > 0.978-R, 



ra s~ 2 , 0-2 X 10~ 2 cm 
2 respectively, with colors ranging from 
Blue contours denote positive values 



lat < 52.5°). Scale ranges are ±9 X 10 
s~ 2 , and 0-8 X 10~ 3 cm s _: 
black/blue to yellow/white 
and green negative. 

components can be represented as follows: 

A ro = Vx (C4>) + Vx , 



(22) 



where £ and x are functions of r and 9. The helioseis- 
mic estimate of Q is obtained from the zonal vorticity 
equation and therefore provides no information about the 
compressible component, \- However, taking the curl of 
equation (l2Tfl) yields an equation for £: 



(VxA m )-0 = -V 2 C + ^ = £ 



(23) 



We solve this equation as described in Appendix C, sub- 
ject to boundary conditions such that the radial acceler- 
ation vanish at the surface (A r = at r = R), and the 
latitudinal acceleration vanish at the base of the NSSL 
(Ag = at r = r s ). The latter condition implies that the 





r > 0.978.R, lat < 52.5° 


r > 0.95.R (extrapolated) 




-8.7xl0~ 8 L 


-2.5xlO" 6 L 




-1.3xlO~ 4 L 


-4.7xl0~ 4 L 




-1.9xl0" 2 L 


-4.2xlO~ 2 L 




8.9xlO" 4 L 


6.7xlO" 3 L 


w t 


-2.0xlO~ 2 L 


-4.9xlO" 2 L 


Ts 


200 days 


590 days 



a Work values are given in terms of the solar luminosity L. 

latitudinal acceleration is relative to the deep CZ. The 
results are shown in Figures [6] and [7] 

To facilitate comparison between the terms, the plot- 
ting range in Figure |5] is limited to that spanned by the 
local inversions. However, the results for A r and Ag ex- 
tend to higher latitudes and deeper levels because they 
are based on the global inversions in Figure [1] This 
greater latitudinal extent is reflected in Figure [3 which 
also serves to illustrate the relative amplitudes and signs 
of each term. 

The amplitudes of Ag and A$ peak at mid-latitudes 
while A r is radially inward at the equator and outward 
above latitudes of about 26°. As noted above with regard 
to Figure 0Jj, the zonal acceleration is negative, tending 
to decelerate the rotation rate in the NSSL. The positive 
sign of Ag indicates an equatorward direction, tending to 
oppose the poleward meridional flow, thereby maintain- 
ing the rotational shear against Coriolis-induced circula- 
tions. The reversing sign of A r , with an inward orienta- 
tion at the equator also tends to oppose the meridional 
flow. This is consistent with the conceptual framework 
put forth in ^and |g] 

The amplitudes of Ag and A^ are comparable, which 
is remarkable because they are approximated by means 
of very different data sets and analysis techniques and 
represent distinct dynamical balances. The amplitude of 
the vertical acceleration is roughly an order of magnitude 
smaller than the horizontal components. 

We can estimate the total rate of work done by each 
of the turbulent stress components as follows: 



W r = 



p s (vi) AidV 



(i = r, 9, 4>) (24) 



where V is the volume of the NSSL and p s is again the 
spherically-symmetric density given by Model S. The re- 
sults are listed in Table 1. 

In the first column of Table 1, the volume V corre- 
sponds to the region spanned by the local helioseismic 
inversions and plotted in Figure [S] In the second col- 
umn we extrapolate this to the entire NSSL, spanning 
0.95i?-1.00i? in radius and pole-to-pole in latitude. The 
poleward extrapolations are based on the analytic fits 
expressed in equations (|16j) and (|17|) , using global inver- 
sions for f2 and local inversions for (vg) . 

When extrapolating A<p in particular, we compute V£ 
in equation (1141) from global inversions interpolated onto 
the local helioseismic grid for r > 0.977i?. We then ex- 
trapolate the expansion coefficients c\(r) and c$(r) that 
appear in (|17|) . as shown in Figure in order to obtain 
(vg) and (v r ) down to 0.95i?. Slopes are chosen to flatten 
out below the matching point at r = 0.977R, anticipating 
an eventual flow reversal deeper in the CZ (we choose C3 
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Fig. 7. — Acceleration terms A r (dashed line, blue in online version), Ag (solid black line), and (dot-dashed line; red in online version) 
as in Figure [6] but here plotted together versus latitude for several different radii, as indicated. The radial component A r is multiplied by 
a factor of 10 for clarity of presentation. Note that we only show the A^ curve out to 52.5° latitude because that is the extent of the local 
inversions used to compute it. 

by means of transporting angular momentum through 
the surface S. In a stationary state the mean mass flux 
is divergenceless, so equation ([2]) implies F = — (pv m ) C. 
The resulting values of W g and W t are listed in Table 1 
and indicate that the latter dominates W^, accounting 
for its negative sign. Thus, most of the work done by 
turbulent stresses to maintain the NSSL against advec- 
tion by the meridional flow involves transporting angular 
momentum out of the layer. 

Our estimate for the angular momentum flux, F, also 
allows us to compute a spin-down time scale for the 
NSSL: 

JyPZdV 

Ts = jF^S ■ (26) 

As indicated in Table 1, this suggest a value of about 200 
days for the region spanned by the local inversions and 
about 590 days for the entire NSSL. This is comparable 
to the ventilation time 
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Fig. 8. — Radial variation of the expansion coefficie nts ci(r) 
(solid line) and cs(r) (dashed line), defined in equation l|17jl . For 
r > 0.977.R (indicated by the vertical dotted line), the curves are 
based on fits to local helioseismic inversions. The extrapolations 
for r < 0.977/? are chosen to have continuous values and derivatives 
but flatter slopes as they extend deeper. 

to flatten out somewhat more rapidly than c\ in order to 
avoid high-latitude counter-cells resulting solely from the 
extrapolation) . The numbers in Table 1 are not expected 
to be greatly sensitive to uncertainties in the meridional 
flow structure, unless the flow reversal actually occurs 
within the NSSL (S3]). 

The values in Table 1 indicate that the zonal force 
(corresponding to the net axial torque J 7 ) does the 
most work, accounting for 2-4% of the solar luminosity. 
The work done by the latitudinal force Ag is roughly two 
orders of magnitude smaller and that done by the radial 
force A r is two to three orders of magnitude smaller still. 

If the net axial torque T can be represented as the 
divergence of an angular momentum flux T = — V-F, 
then the zonal work can be expressed as 



^— - 630 days 



(27) 



W 6 



nF-ds = w a + w t 



(25) 



where S is the bounding surface of V, and dS is directed 
normal to the surface. Here W g is the work done against 
the angular velocity gradient and Wt is the work done 



where Vg ~ 20 m s -1 is a typical meridional flow am- 
plitude. This correspondence is not surprising since we 
estimated the turbulent angular momentum flux F based 
on the meridional flow. 

4.4. Variability and Uncertainty 

As mentioned in Sjl] and as demonstrated in ^4.11 (Fig. 
[3]), estimates of the meridional flow based on different 
observational data sets and different analysis techniques 
can vary substantially. Some of this variability can be 
attributed to measurement uncertainities but some is un- 
doubtedly due to the intrinsic variability of the flow itself. 
In addition to the meridional flow, equations ([2]) and © 
also involve gradients in the angular velocity Q which are 
also subject to measurement uncertainties and intrinsic 
variability. In this section we address what implications 
this variability and uncertainty has with regard to the 
principle arguments, results, and conclusions presented 
in this paper. 
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We begin with the physical processes discussed in sec- 
tions [5] and [3) The predominant dynamical balances we 
advocate for the NSSL, represented by equations © and 
(fTTj) , are expected to hold regardless of the detailed am- 
plitude and structure of the mean flows. They assume 
only that well-defined, persistent zonal and meridional 
flows exist when averaged over time intervals longer than 
a rotation period (of order one month) and that turbulent 
(i.e. Reynolds) stresses dominate over baroclinic forces 
{Q » B) in the NSSL. 

Global rotational inversions indicate that variations in 
£1 over the solar cycle are not more than a few percent 
I Thompson et al.l (|2003f) . Although D, gradients may vary 
more than this (see below), they are not large enough to 
reverse the sign of dC/dX. In order to do so, relative vari- 
ations in the rotation rate Af2/f2 would have to exceed 
2D/R where D is the length scale of the shear. Taking D 
to be the extent of the NSSL ~ 0.05i? implies variations 
with time of the zonal flow of at least several hundred 
m s _1 . Such large temporal variations can be ruled out 
based on photospheric observations and helioseismic ro- 
tational inversions (torsional oscillations are more like 
20 m s _1 ). Likewise, the uncertainty in p (obtained from 
Model S) is constrained by helioseismic structure inver- 
sions and is thus certain to within a few percent. 

Thus, the principle source of uncertainty in the quan- 
titative estimates presented in sections 14.21 and 14.31 is 
the meridional flow. Relative to zonal flow inversions, 
meridional flow inversions are subject to both greater 
measurement errors and greater intrinsic variability, by 
virtue of their weaker amplitude and asymmetric struc- 
ture with respect to the equator (which precludes global 
inversions) . Note that this source of uncertainty only ef- 
fects our estimate for J- and related quantities such as T 
in Fig.HJi, in Figures HJ: and [71 and W^, W g , W t , and 
t s in Table 1. Our estimate for Q and related quantities 
(A r , Ag, W r , Wg) depends only on global rotation inver- 
sions, which are generally more reliable, particularly at 
high latitudes and large depths below the photosphere. 

As noted in §4.21 the radial flow (v r ) makes a negli- 
gible contribution to our estimate of J 7 , even when the 
variability and uncertainty in the meridional inversion is 
taken into account. We have verified that this is consis- 
tent with our data by computing T for each of the years 
shown in Figures H and [JJ (1996-2004). The amplitude 
of the radial component never exceeds 2xl0 7 cm 2 s _1 , 
which is 5% of the saturation value of the color scale 
used in Figure [Hi. The average amplitude is less than 
one percent of this saturation value. 

The time variation of our (vg) inversions varies greatly 
with latitude and depth, as shown in Figure [31 At lati- 
tudes less than 40° the standard deviation of these mea- 
surements is less than 50%. This is also the case for 
our estimate of T , confirming that (vg) is the principle 
source of uncertainty and variation. At latitudes greater 
than 40° and radii below 0.997? the inversions become 
less reliable and the standard deviation exceeds 50%. 

However, we remind the reader that our quantitative 
estimates for T and related quantities in 14.21 and 14.31 
are based on an analytic fit to the 1996 inversions, as 
shown by the dashed line in [3] (1996-2004). This is a 
smooth poleward flow with an amplitude of 18-21 m 
s _1 at all measured radii r > 0.977-R. Although the 



amplitude of this flow is debatable, this general profile 
is consistent with the inpretation put forth by several 
authors who have argued that the spatial and tempo- 
ral variation of the meridional flow can be interpreted 
as a persistent poleward background flow plus a time- 
dependent co mponent associated with c yclic magnetic 
activ i ty (e.g. [S nodgrass & Dailev 1996; Ba su fc Antial 
l20Trl lHathawav fe Rightmirdl201lD . We are interested 
mainly in the background flow which is present through- 
out the solar cycle and is most prominent during solar 
minimum. 

We emphasize that time averages of at least 2-3 months 
are needed to properly assess the mean flows we are con- 
cerned with here. That is the time scale over which 
the dynamical balances in equations © and (jTTJ) are es- 
tablished by means of the Coriolis force. Estimates of 
the meridional flow vary significantly from one Carring- 
ton rotation to the next, independent of the measure- 
ment technique. For example, the results of Hathaway & 
Rightmirc (2010) based on feature tracking indicate that 
the amplitude of the dominant Legendre component can 
change by as much as 5 m s^ 1 over the course of one 
year. If one includes the higher-order components, the 
variation can be even larger, of order 5-10 m s" 1 at the 
higher latitudes. 

The importance of temporal averaging is evident when 
comparing our meridional flow inversions to those of 
Basu & Antia (2010), which are based on the same MDI 
data sets but each averaged over a single Carrington rota- 
tion. The natural variability in our measurements should 
be reduced from those of Basu & Antia by roughly a fac- 
tor of 1.5 just due to our longer averages. Furthermore, 
our estimates are obtained by averaging over both lon- 
gitude and time. At any given instant in time, we make 
189 separate flow determinations scattered across the so- 
lar disk. All measurements made at the same latitude are 
then averaged together (15 determinations at the equator 
and 7 at the highest latitudes). The process is repeated 
on a daily basis and the results further averaged in time. 
Basu & Antia, on the other hand, only make estimates 
along the central meridian. So at high latitudes, even if 
the two teams average over a similar duration, we have 7 
times more data going into the average. This high degree 
of temporal variability is the reason we used a representa- 
tive smooth fit to the helioseismic measurements instead 
of directly using the measurements themselves. 

According to our inversions, the amplitude of 
the background poleward flow is approximately 20 
m s _1 . Other estimates based on various tech- 
niques including heliseismic inversions, Doppler mea- 
surements, and feature tracking generally r ange from 
10-20 m s" 1 (iSnodgrass fc Dailevl Il996t iHathawayl 



1991 iHaberetal.l 120021: iZhao fc Kosovichevl 1 2004 ; 



Hindman et al l l2004MGonzalcz-Hcrnandez'et al.1 2006 



Ulrich 2010; B asu fc Antiall2010l : lHathawav fe Rightmird 
2010L 1201 If ). If the mean meridional flow velocity were 



more like 10 m s -1 as suggest ed by feature tracking (e.g. 
lHathawav fc Rightmirdl2010f L then our quantitative es- 
timates for T in Figure [4ji and A^ in Figures and 
would be a factor of two too large. Furthermore, if 
the actual amplitude of {vg) were to drop to ze r o nea r 
the base of the NSSL as suggested by lHathawav! (|2011f) . 
then our extrapolated values for W^,, W g , and Wt listed 
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in Table 1 would be overestimated. 

With these variations and uncertainties in mind, we 
believe that our quantitative estimate for W$ given in 
£14.31 is reasonable. In particular, since the meridional 
flow speed {v$) is not likely to exceed 20 m s _1 through 
much of the NSSL (r > 0.95i?), then the extrapolated 
value of W<j, is unlikely to exceed four percent of the 
solar luminosity, L. Moreover, since the average value of 
(v e ) throughout the NSSL (r > 0.95R) is not likely to be 
much less than 10 m s _1 , then is unlikely to be less 
than 1-2 percent of L. 

Our estimate for Q is based solely on global rotation in- 
versions from 1996. Still, we can use the local inversions 
for the zonal flow shown in Figure [5] to estimate the tem- 
poral variation. This yields a standard deviation of less 
than 15% for latitudes less than 40°. The maximum stan- 
dard deviation at the highest latitudes measured, 52.5° 
is 34%. This may be regarded as an upper limit to the 
uncertainty in the quantitative estimates for Q, A r , Ag, 
W r , and Wg in sections 14.21 and 14.31 We do not consider 
the time variation of global rotation inversions here but 
we expect they would imply smaller variations of Q by 
virtue of their greater accuracy at higher latitudes and 
lower radii. 

5. THE NATURE OF TURBULENT TRANSPORT 
5.1. Turbulent Diffusion and C Mixing 

Perhaps the simplest paradigm one can envision to po- 
tentially account for the existence of the NSSL is that 
of turbulent diffusion. In particular, one might suppose 
that as the length and time scales of the convection de- 
crease drastically approaching the photosphere from be- 
low, the effects of rotation and spherical geometry be- 
come unimportant, so the convection becomes more ho- 
mogeneous and isotropic, at least in the horizontal di- 
mensions. Furthermore, given the small scale of the con- 
vective motions relative to the scale of the differential ro- 
tation and meridional circulation, one might expect them 
to suppress shear. The decrease in fi near the surface 
may then be attributed to the vanishing of the Coriolis- 
induced velocity correlations that sustain the differen- 
tial rotation deeper in the convection zone (in mean-field 
parlance, this would be the non-diffusive component of 
the Reynolds stress tensor, known as the A-effect; e.g. 
lRiidigerlH989l ). The NSSL would then be maintained 
by the outward diffusion of angular momentum from 
the deep convection zone together with the advection 
of angular momentum by the meridional flow. Is this 
paradigm consistent with observations? 

An alternative paradigm introduced briefly in SJTJ is 
that photospheric convection tends to conserve angu- 
lar momentum locally. This has been found repeat- 
edly in numerical simulations of rotating convection in 
parameter regimes wit h weak rotational influence (high 
Ross by number, see : iGilmanl 1 19771: iGilman fc Foukall 
19791: lHathawavlll982t IDeRosa et al.ll2002t lAurnou et al.l 
2007; Augustson et al. 1 1201 lh . In the absence of ex- 
ternal forcing or boundary influences, turbulent mixing 
tends to establish a rotation profile such that C is uni- 
form throughout the domain in question (17 cx A~ 2 ). 
This is particula r ly evi dent in the recent simulations by 
lAugustson et al.l (|201 lh that employ open boundary con- 
ditions in the vertical, permitting convective motions to 



pass through the boundaries. 

In this section we consider these two paradigms for tur- 
bulent transport and investigate whether they are con- 
sistent with the helioseismic inversions discussed in <2J 
We will focus primarily on the angular momentum trans- 
port, represented by J-, and its implications with regard 
to the fi profile. However, we emphasize again that 
the f2 profiles considered here can only be maintained 
with the help of corresponding meridional forces Q or 
B; these meridional forces are required to maintain rota- 
tional shear dfl/dz regardless of the nature of J- '. 

For both paradigms, we express the net axial torque 
in terms of the divergence of a turbulent angular mo- 
mentum flux: J- = —V-F. This is consistent with the 
hyperbolic nature of the compressible MHD equations 
and reliably captures the various components of J 7 , in- 
cluding the Reynolds stress, Lorentz force, and viscous 
diffusion (Appx. A). 

In the case of turbulent diffusion, the form of the angu- 
lar momentum flux F is analogous to that for molecular 
diffusion but the effective turbulent viscosity v t is many 
orders of magnitude larger than the molecular value. 
This yields 

T = -V-F* = V- {pv t \ 2 Vfi) . (28) 

Note that the direction of the angular momentum flux is 
down the gradient of fi, Ft cx —Vfi, tending to suppress 
rotational shear. By contrast, in our second paradigm, 
turbulence tends to mix angular momentum so we expect 
the flux to be down the gradient of C: 

T = -V-F, = V- {pv a VC) . (29) 

where F a cx — V£. Given the C profile in Figure 16, 
equation (|29|) produces an angular momentum flux that 
is radially inward and poleward, as in high- Rossby n um- 
ber convection simulations (e.g. IDeRosa et aT1l2002| ). 

We emphasize that the motivation for equation (|29p 
is purely phenomenological; this form is suggested by 
simulations of turbulent convection at high Rossby num- 
ber. By contrast, equation (|28p can be derived rigorously 
through the techniques of mean-field theory, assuming 
scale separation and that small-scale motions are homo- 
geneous and isotropic. 

We note also that the paper that introduced the con- 
cept of the near-surface shear layer, Foukal & Jokipii 
(1975; hereafter FJ75), touched on both paradigms. 
They represented the turbulent angular momentum flux 
as a viscous diffusion as in equation (|28|) with a constant 
dynamic viscosity pv t . Substituting this into equation 
@ and assuming a cylindrical rotation profile £1 = fi(A) 
yields 

W^ = -V.F t = ^--^-j . (30) 

It is straightforward to show that this equation is the 
same as equation (1) in FJ75. They then consider the 
advection-dominated limit \v\ » v t and argue that the 
flow will conserve angular momentum, yielding a profile 
C = constant. 

However, it is important to note that FJ75 attributed 
the conservation of angular momentum to the convection, 
not to the mean flow. If C represents the axisymmetric 
zonal flow, then only the axisymmetric meridional flow 
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(v\) (assuming (pv\) sa (p) (ua)) will contribute to the 
advection term on the left-hand-side (lhs) of equation 
(001). 

Thus, there are two ways to interpret the FJ75 re- 
sult within the context of the two paradigms considered 
here. The first interpretation is that the turbulent angu- 
lar momentum flux is diffusive in nature [eq. (|28j) ] and 
the homogcnization of angular momentum C ~ constant 
is achieved by means of the meridional flow. The second 
interpretation is that the (v\) and C = A (v<p} terms on 
the lhs of equation (|30|) represent not mean flows but 
rather convection; that is, replace the mean velocities 
with fluctuating velocities. When averaged over longi- 
tude and time, the net turbulent stress will then tend to 
homogenize £, as expressed in (|29l) . 

Yet, both interpretations are incomplete in the sense 
that they do not take into account the observed pole- 
ward merdional flow which clearly crosses C contours as 
discussed in SJH The meridional flow does not conserve 
angular momentum on its own ((pv,„) -V£ 7^ 0), nor 
does the convection (F 7^ 0). Both must contribute to 
the subtle dynamical balances in the NSSL. 

5.2. Homogeneous Solutions and Stability 

In ij5.ll we suggested two potential idealized forms 
for the turbulent angular momentum flux, expressed in 
equations ([28| and (|29| . We now ask whether these 
paradigms are supported by helioseismic inversions. In 
this section we neglect the meridional flow and consider 
only the homogeneous solutions discussed in H2.2[ equa- 
tion 0. Namely, we calculate the fi profile that would 
result in no net axial torque, J- — and we ask whether 
this profile resembles the f2 profile deduced from helio- 
seismology. Such solutions satisfy our zonal momentum 
equation ^ for (/9v m ) = 0, as they must if they are to 
describe a steady state. We consider the influence of a 
meridional flow in £15.31 

We begin with the case of turbulent diffusion and 
we assume for simplicity that the density-weighted (dy- 
namic) diffusion coefficent pv t is constant. Anisotropic 
and inhomogeneous diffusion coefficients will be consid- 
ered in M5.41 From equation (|2"51) . the condition T = 
then requires 



(31) 



V- (\ 2 Vtt) = (turbulent diffusion) 



Similarly, if we assume pv a is constant, equation 
yields 

V 2 £ = (£ mixing). (32) 

In order to obtain solutions to equations (pTTj) and (|32p 
we must specify boundary conditions. Thus, to proceed, 
we must keep in mind the context in which these equa- 
tions are proposed to be valid. They are intended to rep- 
resent angular momentum transport in the NSSL alone; 
angular momentum in the deep convection zone (CZ) 
must be very different in order to sustain the solar differ- 
ential rotation. Thus, to take into account the coupling 
between the NSSL and the deep CZ, we specify boundary 
conditions at the base of the NSSL. 

It is often stated that the fi contours at mid latitudes 
in the bulk of the CZ are nearly radial. Close scrutiny of 
Figure Q}z reveals that f2 contours are not strictly radial; 
rather, they are tilted slightly toward the rotation axis so 



dil/dr > in the bulk of the convection zone. However, 
in the NSSL, dfl/dr < 0. Thus, we define the base of 
the NSSL, r s , as the radius at which the spherically- 
averaged radial 51 gradient passes through zero. This 
yields r s = 0.946. 

Thus, we solve equations (|3"Tj) and (|32| in the region 
r s < r < i?, subject to the boundary conditions dfl/dr — 
and O = f2 s (8) at r — r s . Here f2 s (8) is an analytic fit 
as in equation (1161) to the helioseismic inversions shown 
in Figure [T] interpolated to r — r s . For further details 
on the boundary conditions and for the analytic solution 
of equations (f3"T|) and f[32]) see Appendix D. The results 
are plotted in Figure |9 

It is immediately apparent in Figure |H] that the ac- 
tual il profile in the Solar NSSL inferred from helioseis- 
mology is steeper than suggested by either of the sim- 
ple paradigms considered here. Not surprisingly, viscous 
diffusion tends to suppress shear so the equilibrium 
profile in this case is nearly independent of radius. How- 
ever, the latitudinal differential rotation is still promi- 
nent, roughly the same as at r = r s as a consequence 
of the thinness of the layer (coupled with the boundary 
conditions). Furthermore, a weak negative radial shear 
(d£l/dr) is established at high latitudes. This is a geo- 
metric effect associated with the inefficiency of viscous 
transport near the rotation axis (A = 0). 

The profile corresponding to C mixing [eq. (|32[) ] is 
steeper than the diffusive profile, again as expected, since 
the direction of the flux (cx — V£) has a component that 
is radially inward. However, note that, within the con- 
text of the NSSL, the concept of turbulent convection 
mixing angular momentum is not the same as turbulent 
convection establishing C = constant. The thin-shell ge- 
ometry and the coupling of the NSSL to the rotation 
profile in the deep CZ (modeled here by means of our 
boundary conditions) preclude a complete homogeniza- 
tion of C even if the local convective transport were to 
exhibit that tendency. 

At high latitudes, the f2 profile implied by C mixing 
is steeper and is more comparable to the helioseismic in- 
versions. Both are less steep than the profile that would 
arise from a complete homogenization of C. However, 
the helioseismic profile does approach homogenization in 
radius within about 7 Mm of the photosphere. In other 
words, the upper bound on the magnitude of dfl/dr near 
the photospher e appears to b e set by the Rayleigh sta- 
bility criterion (|Tassoulll 19781) 



dC 
dX 



> 



(V5 = 0) 



(33) 



Note that equation (|33l is only strictly valid under the 
assumption that the stratification in the convection zone 
is approximately adiabatic (V5 = 0). This is of course 
an oversimplication since the stratification in the NSSL is 
thought to be substantially superadiabatic. A more rig- 
orous analysis indicates that the Rayleigh criterion (|33l) 
is part of a more general formulation of the Schwarzchild 
criterion for convective stability, which is clearly violated 
in a convection zone, essentially by definition. An inde- 
pendent stability criterion requires cos^VSx V£) > 0. 
However, if S isosurfaces are predominantly horizontal, 
then this provides no constraints on the radial shear in 
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Fig. 9. — Rotation rate profiles are shown versus radius, spanning the upper convection zone (0.92i? < r < Ft) for latitudes of (a) 0°, (6) 
20°, (c) 45°, and (d) 69°. Dot-dashed lines (red in online color version of Figur e) a nd da shed lines (blue in online version) represent viscous 
diffusion and C mixing respectively, obtained by solving equations equations H3U and (132 I t subject to the boundary conditions discussed 
in the text. Solid lines (black in online version) represent the helioseismic inversions shown in Fig. [T] and vertical dotted lines indicate 
the matching layer, r = r s . The sloped dotted line (green dashed line in online version) represents an angular momentum profile that is 
independent of depth: C = Cr(8) where Cr(0) is the surface value (at r = Ft). 



the NSSlB, dVt/dr. Both criteria arise from the Solberg- 
H0iland stabil ity analysis as discussed, for example by 
iTassoull (|1978|) . By considering the implications of equa- 
tion (1331) here, we are essentially separating out inertial 
effects from thermal effects, which may operate on differ- 
ent time scales. The potential relevance of this separa- 
tion to the NSSL is supported by the remarkable corre- 
spondence between the maximum radial gradient of the 
helioseismic f2 profile and the slope implied by equation 
(|33| . as illustrated in Figure [9] 

Thus, with these caveats, equation (|3"31 implies that 
the radial gradient near the photosphere is what it is 
because anything steeper would be unstable. Although 
this is indeed a compelling argument, it is well known 
that in the presence of a weak magnetic field the Rayleigh 
stability criterion is replaced by the stability criterion as- 
sociated with the magneto- rotational instability (MRI). 
Thus, equa tion (1551) is replace by the condition that 
dfl/dX > (|BalbusH1995f) . This condition is clearly vio- 
lated in the NSSL. 

Why might the rotation profile be limited by the hy- 
drodynamic Rayleigh criterion (|33[) yet violate the MRI 
stability criterion? There are two potential answers to 
this question. The first is that the MRI stability analysis 
may n ot be applicable in the NSSL. According to[Balbus 
(1995), the assumptions that underlie the MRI crite- 
ria are valid for field strengths in the range AitpxQ « 
B 2 « 4irpR 2 il 2 , where \ is the transport coefficient (in 
units of cm 2 s^ 1 ) corresponding to the dominant diffu- 



sive process. Dissipation in the NSSL is likely dominated 
by radiation and ohmic diffusion. The detailed physics 
is complicated and very sensitive to depth but never- 
theless, rough estimates for the corresponding diffusion 
coefficients give k ~10 4 -10 5 cm 2 s _1 and r\ ~ 10 6 -10 7 
cm 2 s _1 . The k value is derived from sola r structure 
Model S (jChristensen-D alsgaar d et al.lll996D and the r\ 
value follows from the Spitzer expression for a fully ion- 
ized Hydrogen plasma 77 - 8 x 10 13 T" 3 / 2 (jSpitzerlll962t) . 
Using p ~ 10~ 3 g cm -1 , we estimate that the MRI anal- 
ysis of iBal bus (1995D is valid for field strengths ranging 
from less than 1 G to more than 10 4 G. The typical field 
strength in the NSSL is likely to lie within these bounds. 

The second possibility is that the hydrodynamic 
Rayleigh instability is somehow more robust or more ef- 
ficient than the MRI. Thus, the NSSL may indeed be 
unstable to MRI but the time scale of the instability (of 
order the rotation period ~ 28 days) is longer than the 
time scale over which the shear is established by tur- 
bulent stresses (of order the convective turnover time ~ 
5 min - 1 day). Although this is plausible, the same 
argument would also in principle apply to the Rayliegh 
criterion so it is unclear why equation (|33[) must be satis- 
fied while the MRI criterion is not. The depth over which 
the marginal slope is achieved, roughly within ~ 7 Mm 
of the photosphere, suggests that granulation may play a 
role. In any case, this is an interesting issue that should 
be explored further with the help of MHD convection 
simulations. 



In fact, it simply reflects convective instability since 
cos8(dS/dr)(dC/de) < 



5.3. Balancing the Meridional Flow 
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Fig. 10. — Differential operators applied to the (smoothed) solar 
f2 profile inferred from helioseismology (Fig. [T}, corresponding to 
(a) viscous diffusion V-(A 2 V!2.) and (b) C mixing, V 2 (A 2 Q«). 
Solid lines denote negative values, dotted lines denote zero and 
positive values. The scaling is arbitrary since the corrcspoding 
torques depend on the unknown values of vt and v a . 



The analysis presented in §5.21 assumes that the net 
torque F is zero, which is clearly not the demon- 
strated in 21 observations imply that F is negative. 
Thus, we must take this into account if we are to prop- 
erly assess whether our two simple paradigms for angular 
momentum transport can adequately account for the ob- 
served rotation and meridional flow profiles. 

Stated another way, observations imply that the merid- 
ional flow is supplying angular momentum to the NSSL, 
tending to speed up the local rotation rate, while the tur- 
bulent transport must be removing angular momentum 
in order to maintain a stationary state. In order to assess 
whether turbulent diffusion or £ mixing can provide the 
requisite transport, we must compute F from equations 
(|2"8|) and (j2T)]) based on the rotation profile inferred from 
helioseismology £7, (r, 9) and ask whether it matches the 
inferred net axial torque shown in Figure |4]2. 

Of course, we cannot compute F explicitly because we 
do not know a priori the transport coefficients v t and v a . 
However, if we assume as in £15.21 that pv t and pv a are 
constant, then the corresponding torques, F% and F a will 
be proportional to the appropriate differential operators 
applied to the solar rotation profile, namely 



F t oc V (A 2 Vft») and F a oc V 2 (A 2 ft» 



(34) 



The right-hand-side of each of these equations is plotted 
in Figure [TOl 

The immediate impression from Figure [10] is that the 
two results look nearly identical. This is a consequence 
of the steepness of the radial il gradient and the thin- 
shell geometry. Thus, the dominant contribution to both 
operators is the second radial derivative: 



v (a 2 v^*) « v 2 (x 2 n, 



dr 2 



(35) 



As shown in Figure H] the solar ft profile is nearly flat 
at the base of the NSSL (dQ/dr ~ 0) and steepens with 



increasing radius, reaching its maximum slope near the 
photosphere. Thus, d 2 ft/dr 2 < and the quantities 
shown in Figure [TU] are predominantly negative in the 
NSSL (r > OMR). 

The negative sign of the quantities plotted in Figure [TU] 
bodes well for the viability of our two simple paradigms 
for angular momentum transport. Both viscous diffusion 
and £ mixing would tend to flatten out the negative cur- 
vature of the solar f2 profile (d 2 ft/dr 2 < 0). This would 
tend to decelerate the local rotation rate (F < 0) and 
may thus serve to remove the angular momentum sup- 
plied to the NSSL by the meridional flow. The efficiency 
of the transport scales with the coefficients v t and v a and 
could be calibrated to give the proper net torque, at least 
in an integrated sense. 

However, the profiles shown in Figure [10] are clearly 
different than the net axial torque profile inferred from 
helioseismic measurements, shown in Figure HJi. In par- 
ticular, the idealized profiles in Figure [TU] peak at the 
equator while the helioseismic profile in Figure 0]z peaks 
at mid latitudes, with a weak positive signal at the equa- 
tor. Thus, our simple paradigms appear to be inconsis- 
tent with helioseismic inversions. However, it may plau- 
sibly be argued that the assumption of homogeneous, 
isotropic coeffients pv t and pv a is unrealistic and unnec- 
essary. In the next section we consider whether this can 
redeem our simple paradigms. 

Before proceeding, we note in passing that the C mix- 
ing operator shown in Figure I10t > suggests a reversal in 
the sign of the torque F at high latitudes. Since the lo- 
cal helioseismic inversions only extend up to latitudes of 
about 50°, we cannot determine whether a similar sign 
change occurs in the Sun. However, some measurements 
of the meridional flow suggest that there may be per- 
sistent high-latitude counter-cells, which are particularly 
apparent at solar minimum when magnetic eff ects can 
be more easily separated out (e.g. lUlrichl 120101 see also 
Fig. [3]). Under the justified assumption that the angu- 
lar momentum continues to decrease toward the poles 
(dC/d9 < in the northern hemisphere), equation ([2]) 
implies that a high-latitude counter-cell would indeed 
correspond to a change in sign of the net axial torque 
F (such that F > 0). 

5.4. Anisotropy and Inhomogeneity 

The idealized rotation profiles and net axial torques 
considered in ^5.21 and i j5.3l are only valid for constant 
density- weighted transport coefficients pv t and pv a . If 
one allows for inhomogeneous, anisotropic diffusion ten- 
sors then can the two paradigms considered here provide 
a better fit to helioseismic inversions? 

This is to some extent a tautology; one might expect 
that one can construct a diffusion tensor to reproduce 
any arbitrary f2 and F profiles so the results of this anal- 
ysis would not be very enlightening. However, we will 
demonstrate that turbulent diffusion can be ruled out as 
a viable paradigm even if it is anisotropic and inhomo- 
geneous. Mixing angular momentum, on the other hand, 
may be more promising. 

Note that, with regard to anisotropy, we are referring 
explicitly to the off-diagonal components of the viscous 
stress tensor u t or its C- mixing analogue, v a . This should 
not be confused with the off-diagonal components of the 
Reynolds stress tensor, which are nonzero even for con- 
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stant scalar values of v t and v a . Indeed, the off-diagonal 
components of the Reynolds stress tensor are essential in 
order to account for the maintenance of mean flows. 

Perhaps the most conspicuous fault with the profiles 
derived in £|5.2l is that they admit an angular momentum 
flux through the solar surface. In the Sun, the sharp drop 
in density near the surface precludes any significant an- 
gular momentum flux through the photosphere on a time 
scale comparable to the convection turnover time or the 
rotation period. Angular momentum is lost through the 
solar wind but this loss rate is many orders of magni- 
tude smaller than the rate at which angular momentum 
is continually circulated through the NSSL by convection 
and meridional flows. 

Thus, we consider equation (|28|) again but we now al- 
low for a diffusion coefficient that depends on radius. 
Furthermore, we allow for anisotropic transport so we 
can regard v t as a tensor. For the time being we will 
neglect the off-diagonal elements of the tensor so we can 
express v t in terms of vertical and horizontal component 
v tv and v th . 

For pv tv — constant we found in i)5.2l that the resulting 
profile was too shallow to account for the helioseismic 
fl profile. Can diffusion give rise to a steeper profile? 
The answer is yes, but only if the vertical diffusion pv tv 
increases with radius {d(pv tv ) / dr > 0). This would pro- 
vide an outward angular momentum flux that becomes 
larger closer to the solar surface, providing the requi- 
site flux divergence to decelerate the NSSL. However, as 
pointed out in the previous paragraph, this cannot be 
sustained all the way to the photosphere where F and 
thus pv v t must drop to zero. 

To demonstrate that even an anisotropic, inhomoge- 
neous diffusion cannot account for the NSSL consider 
a closed volume V bounded from below by the surface 
r = rf, and from above by the photosphere R. Latitu- 
dinal boundaries are at the equator and at an arbitrary 
colatitude 9q 1 located in the northern hemisphere close 
enough to the equator to be accessible to helioseismic 
inversions, say 7r/6 < 9q < tt/2 (corresponding to a lati- 
tude between and 60°). We now integrate the net axial 
torque T over the volume, allowing for anisotropic trans- 
port and assuming F = at r = R. The symmetry of the 
Q profile implies no angular momentum flux through the 
equatorial plane {dVL/39 = 0) so the integrated torque is 
given by 



/ TdV = - [ F-dS 
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where S is the bounding surface, and the first and second 
terms on the right-hand-side are to be evaluated at r = rb 
and 9 = 9o respectively. 

As discussed in Sj4] and £)5.3[ the integral in equation 
(|36|) must be negative. However, if we choose r& to lie 
within the NSSL, then dil/dr < and the first term is 
positive. In other words, angular momentum flux into 
the NSSL through the lower surface r — n must be 
shunted poleward in order to be consistent with a nega- 
tive (or zero) net torque J v TdV . Although the latitu- 
dinal flux is indeed poleward (dQ/89 < at 6 = 9b), it 



cannot be efficient enough to maintain a negative 8Q /dr 
throughout the NSSL. As n, approaches R from below, 
the second term in brackets can be approximated by 
pvhR 2 (R — rb)(dil/d9). In order to transport the req- 
uisite flux, the viscosity anisotropy Vh/v v would have to 
increase radially without bound, becoming infinite at the 
photosphere. 

This argument can be readily generalized to rule out 
any arbitrary outward flux F-f > 0. The implication, 
then, is that the angular momentum flux in the NSSL 
must be radially inward. One can in principle salvage 
the diffusive paradigm if one invokes negative diffusion 
or the off-diagonal elements of the turbulent viscosity 
tensor, namely an inward angular momentum flux that 
is proportional to the latitudinal shear. However, such 
prescriptions seem rather contrived. 

Thus, we can confidently say that the turbulent trans- 
port in the NSSL cannot be adequately modeled as a 
turbulent diffusion. How does our other paradigm fare, 
namely that of mixing angular momentum? 

We can see immediately that this paradigm is more 
plausible because the angular momentum transport is 
indeed inward F-f < 0. A steeper dfl/dr profile could 
then be achieved by means of a diffusion coefficient pv a 
that decreases with radius (pv a )i approaching zero at 
the photosphere, as required by the condition F = 
at r = R. Thus, the £-mixing paradigm is consistent 
with both an inward £1 gradient (dtt/dr < 0) and no net 
angular momentum flux through the photosphere. 

However, in order to be consistent with the form of 
the net axial torque inferred from helioseismology, the 
density-weighted mixing coefficient p = pv a would have 
to have a very particular form. As shown in Figure [Hz, 
the net axial torque is nearly zero at the equator whereas 
V 2 £ has a prominent peak at the equator (Fig. [TUb). In 
order to see how this may be remedied, assume that p is 
isotropic but varies with radius with a scale height much 
less than the solar radius R. Furthermore, assume as in 
equation Q that the C profile is nearly cylindrical. Then 
equation pt))) yields 



T ps p\7 2 £ 



„ . n dpdL „d 2 p 



S ■ (37) 



li p = pv a decreases with increasing radius as sug- 
gested in the previous paragraph, then inserting the he- 
lioseismic value for C = A 2 f2* implies that the first two 
terms are negative at the equator. In other words, the 
sharp decrease of p exacerbates the prominent negative 
net torque T seen in Figure [T7Jb . One way out of this 
dilemma is if the slope flattens out near the photosphere 
so d 2 p/dr 2 > 0. Then the last term in equation (|3T[) 
is positive and correspondence with Figure [4b: becomes 
possible. 

In summary helioseismology sets fairly strict con- 
straints on the nature of the turbulent angular momen- 
tum flux F in the NSSL, including: (1) F must be radially 
inward (F-f < 0), (2) there must be no flux through the 
photosphere (F-f = at r = R), (3) F must diverge 
at mid-latitudes (F < 0) in order to balance the advec- 
tion of angular momentum by the meridional flow, and 
(4) turbulent transport at low latitudes must redistribute 
angular momentum in such a way as to support a steep 
radial Q gradient while minimizing the net axial torque 
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(F = V-F = 0). The turbulent diffusion paradigm, 
equation f2"5|). does not meet these constraints and is 
therefore not a valid model of turbulent transport in the 
NSSL. The alternative ansatz of C mixing, as expressed 
in equation ()29l) , is at least feasible, but it would require 
a fine-tuned, inhomogeneous and/or anisotropic mixing 
coefficient fi — pv a . For example, it could in principle 
be achieved with a vertical mixing coefficient that de- 
creases with radius but flattens out near the photosphere 
(dp/dr < 0, d 2 p/dr 2 > 0). A power law dependence 
(i a r™ with n < may satisfiy such a requirement. 
However, it is more likely that the turbulent Reynolds 
stress in the NSSL is more complex than either of these 
crude, local models. 

5.5. Meridional Momentum Transport 

Despite the conclusion of section §5.41 regarding the 
non-diffusive nature of J 7 , it is reasonable to expect that 
Q may operate essentially as a turbulent diffusion. In 
order to appreciate why this may be the case, consider 
again the time-dependent thought experiment discussed 
in §321 Here we begin with a spherical volume V in uni- 
form rotation and a retrograde torque J- is introduced 
in the surface layers (r > 0.95). This will establish a 
poleward flow (vg) which will steadily increase in ampli- 
tude, striving to establish a cylindrical rotation profile as 
discussed in £13.21 Given the thin-shell geometry of the 
NSSL, one may expect strong vertical gradients dvg/dr 
to be established. Furthermore, given the small scale 
of photospheric convection relative to the mean flow, 
one may expect turbulent stresses to resist such shear- 
ing motions. This would imply that the resulting axial 
differential rotation profile dfl/dz is determined by how 
efficiently small-scale convective motions can mix latitu- 
dinal momentum, (vg), or equivalently, suppress zonal 
vorticity, (w^). 

Thus, under this scenario we may expect that Q will 
be diffusive (down-gradient) in nature, and furthermore, 
that the strain rate tensor will be dominated by vertical 
gradients in the poleward flow. The turbulent transport 
would then be given by 



dr 2 



d 3 



dr 3 



(38) 



where vt is again the turbulent (kinematic) viscosity. 
Note that this expression assumes that the scale of vari- 
ation of the dynamic viscosity \d\a.(pvt) / '<ir| _1 is larger 
than that of the shear |dln |i)g|/cir| _1 . 

Equation (|38j) can in principle provide the momentum 
transport required to balance the Coriolis force associ- 
ated with the axial shear dVl/dz. However, this is con- 
tingent on the amplitude of vg being strongly peaked 
near the photosphere, such that d 3 \vg\/dr 3 < 0. This 
does not appear to be supported by local helioseismic in- 
versions (§31), but the sensitivity and resolution of the in- 
versions is likely not adequate enough to provide reliable 
estimates for third-ord er derivatives of vg (e.g. iBeckersI 
120071: lHathawavll201lD . An alternative is that the dy- 
namic turbulent viscosity pvt decrease with radius while 
the velocity amplitude increase, consistent with no tan- 
gential stress at the photosphere. 

In summary, we can say that if the turbulent transport 
Q does indeed dominate over baroclinic forcing B as sug- 



gested in ^3.3l then it must resist the poleward meridional 
flow induced by T through the Coriolis force. Thus, Q 
must be negative in the northern hemisphere and posi- 
tive in the southern hemisphere. A turbulent diffusion 
or similar mixing process may be adequate, as would a 
more general formulation. 

6. SUMMARY AND CONCLUSION 
6.1. Maintenance of the Solar NSSL 

We have demonstrated that the turbulent angular mo- 
mentum transport in the solar Near-Surface Shear Layer 
(NSSL) is responsible for the persistent poleward merid- 
ional flow but it does not uniquely determine the mid- 
latitude n profile. Rather, the axial rotation gradient 
dfl/dz must be maintained by turbulent stresses in the 
meridional plane. More specifically, a retrograde zonal 
force (axial torque) T < establishes the NSSL and 
regulates the poleward meridional flow while meridional 
forcing regulates the mid-latitude Q profile. Further- 
more, we argue that a transition in the meridional force 
balance from baroclinic to turbulent stresses [B to Q in 
equation (j8|)] may determine the base of the NSSL. 

More generally, we have demonstrated that there is 
close dynamical relationship between the differential ro- 
tation and the meridional circulation and that the struc- 
ture of the NSSL as inferred from helioseismology relies 
on a delicate nonlinear, nonlocal interplay between the 
two. The same physical mechanism (F < 0) that estab- 
lishes negative radial shear dil/dr < also establishes 
poleward flow ((vg) < in the northern hemisphere). As 
suggested by previous authors, we attribute this physical 
mechanism to Reynolds and possibly Maxwell stresses as- 
sociated with the relatively small-scale convection (gran- 
ulation to supergranulation) that permeates the NSSL. 

Throughout our analysis, we have considered the iner- 
tia of the mean flow explicitly, incorporating other phys- 
ical processes in the generalized zonal and meridional 
forcing terms J- and Q which we refer to as turbulent 
stresses. This is intended to clarify the essential physics 
of how the differential rotation and meridional circula- 
tion are coupled in their response to zonal and meridional 
forcing, independent of the detailed nature of this forc- 
ing, which is unknown. The reader may wish to regard T 
and Q simply as the the convective Reynolds stress, since 
this is likely to be their dominant component (Appendix 
A). However, as noted in the introduction, other forcing 
such as large-scale Lorentz forces and viscous diffusion 
can induce mean flows in an analogous way. 

We have demonstrated that turbulent transport in the 
NSSL is non-diffusive in nature and must be directed ra- 
dially inward (§3]). Inspired by numerical simulations of 
turbulent convection at large Rossby numbers, we have 
considered an alternative paradigm for turbulent trans- 
port based on the mixing of specific angular momentum, 
C We have shown that the conservation of angular mo- 
mentum alone cannot account for the existence of the 
NSSL (5j3j) but the form of the rotation profile may be 
consistent with C mixing if the mixing coefficient is in- 
homogeneous and/or anisotropic (§3]). Furthermore, we 
have shown that C mixing is a more general concept 
than simply £-homogenization (constant C) when one 
takes into account the coupling between the NSSL and 
the deep convection zone. Yet, the failure of simple tur- 
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bulcnt diffusion or £-mixing prescriptions demonstrates 
that the NSSL is not simply a passive response to deeper 
forcing; it must be actively maintained by anisotropic 
and inhomogeneous turbulent transport. 

Estimates based on local and global helioseismology 
indicate that it takes 2-4% of the solar luminosity to 
maintain the NSSL against the inertia of the mean flow 
fi j4.3[) . Most of this work is associated with transport- 
ing angular momentum out of the layer (Wt in Table 
1) in order to balance the convergence of angular mo- 
mentum flux into the layer by meridional flow advection. 
The estimated amplitudes of turbulent transport in the 
latitude and longitude directions (Ag and A$) are re- 
markably similar, given the very different way in which 
these two quantities were obtained (the former follows 
from "uncurling" global f2 inversions while the latter in- 
volves estimates of zonal and meridional flows from local 
helioseismology; see £|4.3p . The vertical transport in the 
middle of the NSSL is about an order of magnitude less 
(at the base of the NSSL, r = r s , Ag is zero by assump- 
tion and our local helioseismic inversions provide no in- 
formation on Aff,). The sense of all these terms is such 
that turbulent transport is decelerating the rotation rate 
in the NSSL and opposing the meridional flow. 

Estimates of the spin-down time scale indicate that it 
is similar to the ventilation time scale of about 600 days 
f fc|4.3[) . This is long compared to the turnover time scale 
of convective motions, < 1 day, implying that the net 
turbulent angular momentum transport is rather ineffi- 
cient. 

Finally, we note that the upper limit to the slope of the 
radial 51 gradient appears to be set by the Rayleigh cri- 
terion, equation (1331) . However, it is unclear why this in- 
trinsically hydrodynamic condition should be satisfied in 
the NSSL while its MHD analogue, the stability criterion 
for the magnetorotational instability (MRI) dil/dX > 0, 
is clearly violated. 

6.2. Implications for Numerical Models 

Several authors have sought to investigate the dy- 
namics of the NSSL through numerical simulations of 
solar convection in thin spherical shells or spherical 
segments, placing the lower boundary of the simula- 
tion domain wit hin the upper convection zone, typically 
above 0.9-R (e . g.lGilman k, Foukal|[T979l:lHathawavl ll982l: 
iDeRosa et al"1l2002t lAugustson et all 1201 lft . Relative to 
global simulations spanning the entire convection zone, 
these have the great advantage that smaller scales can 
be resolved so the turbulent transport can be more reli- 
ably captured. Although these simulations have provided 
insight into the dynamics of the NSSL, none have accu- 
rately reproduced the angular velocity profile throughout 
the NSSL. As we have demonstrated here, the dynamics 
of the NSSL involves a delicate balance between small- 
scale turbulent transport, large-scale mean flows, and 
coupling to the deep convection zone (CZ). Even if thin- 
shell models properly capture the turbulent transport, 
they must also capture or otherwise mimic the coupling 
to the CZ if they are to achieve solar-like mean flows. 

A straightforward and common strategy in thin-shell 
models is to impose a solar-like latit udinal differential 
rotat i on on the lower bou ndary (e.g. iGilman fc FoukaH 
119791: IDeRosa et ail 12002ft . If the lower boundary is 
also impermeable, as is often the case, then it is clear 



that the system cannot sustain a poleward flow through- 
out the layer. What implications might this have for 
the differential rotation? As is demonstrated in SJU the 
meridional circulation supplies angular momentum to the 
NSSL while turbulent stresses must remove this angular 
momentum. This balance, together with the meridional 
stresses Q (or B) determine the mean flow profiles. If the 
poleward flow is artificially suppressed by an impenetra- 
ble boundary condition, then there must be alternative 
source of angular momentum to balance turbulent trans- 
port. Without this source, the f2 profile will be adversely 
affected in addition to the meridional flow profile. 

The rotational coupling between the convection zone 
and the NSSL in this class of numerical models is 
achieved by means of viscous diffusion. This is artifi- 
cial in the sense that the viscosity used is many orders of 
magnitude larger than that of the solar plasma. However, 
can this effectively mimic the coupling that is expected to 
occur in the Sun? On the positive side, viscous coupling 
can indeed serve as an angular momentum source, allow- 
ing the integrated turbulent stresses J„ TdV throughout 
the NSSL to be negative, as in the Sun. However, the lat- 
itudinal distribution of the viscous torques is likely to be 
very different than for gyroscopic pumping. If dQ jdr < 
at nearly all latitudes in the solar NSSL as implied by he- 
lioseismic inversions (although these are uncertain pole- 
ward of 70°), then viscous coupling would imply outward 
angular momentum transport everywhere. By contrast, 
the meridional flow would impart angular momentum at 
low latitudes and remove it at high latitudes in such a 
way that the net convergence into the layer is positive. 
This implies a high-latitude angular momentum flux that 
may be up the gradient of f2. One would expect the re- 
sulting 17 profile to be very different than for viscous 
coupling. Mean flow profiles would be similarly sensitive 
to impcrmablc boundary conditions in latitude. 

Thus, in order to properly model the NSSL, numerical 
models must either be deep enough to capture the closed 
meridional circulation, including the poleward flow at the 
surface and the return equatorward flow, or they must 
impose boundary conditions that are conducive to estab- 
lishing solar-like mean flows. These may include open 
boundaries on which the rotation profile CI is specified 
in addition to an imposed net axial torque J 7 , chosen 
to induce a commensurate meridional flow through the 
boundary by means of gyroscopic pumping, as expressed 
by equation J2]). 

Even if the simulation domain is in principle large 
enough to capture the complete, closed meridional cir- 
culation, the delicate balance expressed in equation ([2]) 
implies that the meridional flow is particularly sensitive 
to artificial viscous dissipation. If viscous angular mo- 
mentum transport largely balances the angular momen- 
tum transport by the convective Reynolds stress as in 
many numerical models, then T will nearly vanish and 
the meridional f low may be very different than what oc- 
curs in the Sun (jMiesch et al.ll2008l [2011ft . Regardless of 
the boundary conditions, realistic meridional flow pro- 
files require minimal viscous dissipation. 
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APPENDIX 



APPENDIX A: EXPLICIT EXPRESSIONS FOR T AND Q 

In this Appendix we explicitly identify what is included in the turbulent stress terms J- and Q that are introduced 
in sections 1 2 . II and [3.11 and that are used throughout the paper. 

We begin with the equation that expresses the conservation of momentum in a compressible, electrically conducting 
fluid under the magnetohydrodynamic (MHD) approximation 

P?r + 9 ( v ' v ) v = - VP + P9 + t- (V xB) xB + V-T) (Al) 
at 47r 



Traditional notation is used: v is the bulk velocity of the fluid, p is the density, P is the pressure g = 



—gf is the 
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gravitational acceleration, B is the magnetic field, and T> is the viscous stress tensor, with elements 

1 — val 



T>ij = -2pv 



%i - 3 (V-v) 



(A2) 



where v is the kinematic viscosity. We use spherical polar coordinates (r, 9, <f>) and throughout the bulk of the paper, 
we consider an inertial reference frame. This lets us more gracefully incorporate the differential and uniform rotation 
components into a single, non-uniform rotation profile, Q. However, for the benefit of readers, in this Appendix we 
wish to illustrate explicitly where the Coriolis force enters into this analysis. 

Thus, we can convert equation (|Alj) into a rota ting coordinate system (r, 9, <f>') by writing v = u + QoXcp and 
4> r = <fi + flat. Substituting these changes into (|Alj) yields 



dw 1 
p— + p(u-V) u = —VP + pg- 2pJ7 xu - pfl x (rj xA) + — (VxB) xB + V-T> r (A3) 

Ot 47T 

where f2o = Qoz, A = AA. All derivatives are with respect to the rotating coordinate system so <p may be formally 
replaced by <f> r in the V operators. However, this is not necessary since d/dcf) r — d/d<f> for fixed r, 9, and t. The stress 
tensor XV is the same as T> with v rep lace d by u (and (j> by <fi r ). 

Multiplying the zonal component of (|A3|) by A and averaging over longitude (<j) r ) yields 

^ (pXu^) + (pu m ) -V£ = T = V- [(pXu%) - (XBB^) - puX 2 Vil] (A4) 

where C = X {(u^) + Af2o) = A (ty,) = A 2 0. The right-hand-side is defined as the net axial torque J- '. Thus, it includes 
the Reynolds stress (first term), the Lorentz force (second term), and the viscous diffusion (third term). The Lorentz 

force may be decomposed into a contirbution from mean fields (B) (B^) and a Maxwell stress, /b'P^Y Meanwhile, 

the left- hand- side of (|A4[) includes the Coriolis force and the nonlinear advection (inertia) associated with the mean 
flows, ((v) -V) (v). 

We emphasize that equation (|A4p follows directly from equation (1A3[) with no additional assumptions [although we 
have used the mass continuity equation dp/dt — — V*(pu) in the derivation]. Equation (|A3j) in turn follows dir ectl y 
from (jAip . The only assumption in any of this derivation is the MHD approximation that underlies equation (|A1[) . 
If we make the justified assumption that (pu^) sa (p) {u^) (as, for example, in the anelastic approximation), then the 
first term on the left-hand side of (IA4[) is just d((p) C)/dt. Furthermore, note that u m = v m . We thus obtain equation 

©• m 

Now consider the meridional components of (| A3|) . As is well known, the gravitational and centrifugal terms can 

be expresed in terms of a gradient g + f2oX(f2oXA) = V( v f , s + A 2 57q/2) where ^ g is the gravitational potential. 
Furthermore, we may combine the advection and Coriolis terms as follows 

(u-V) u + 2f2 Xu = uxu + V (-^j ' ( A5 ) 

where u> — Vxu + 2 17q = Vxv is the vorticity relative to the inertial frame, also referred to as the absolute vorticity. 
We can then divide (|A3|) by p, average over longitude, and compute the zonal component of the curl to obtain 

gW) X — = B I g- V{P)xV{p) | [ VPXVP V(^)xV(p) \ 

at dz {p f \ p 2 {p f I 

+ {Vx [(Vx (u m )) x (u m )}}-4> 

+ jvx |(Vxu') xu' + ^-(VxB) xB + p- x V-I> r \|-^ . (A6) 

Again, this equation follows directly from equation (|A3|) with no further assumptions. The second term on the left- 
hand-side, involving fl 2 includes the Coriolis force and the inertia associated with the differential rotation, (u^). The 
term on the right- hand-side involving the mean meridional circulation, (u m ) is estimated to be about two orders 
of magnitude smaller than this. Since we wish to focus on the primary components that contribute to the force 
balance and the inertia of the mean flow, we include the meridional circulation term with the turbulent stress term Q . 
Alternatively, since the kinetic energy density of the convection is at least two orders of magnitude larger than that in 
the meridional circulation, it is justified to neglect the meridional flow term in Q altogether, relative to the Reynolds 
stress. 

The first term on the right-hand-side of equation (1A6[) is the baroclinic term associated with the mean stratification, 
(p) and (P), which we define as B in equation The remaining terms on the right-hand-side of equation (IA6[) 
define the turbulent stress Q . Thus, Q includes residual baroclinic forcing involving fluctuating density and pressure 
components p' and P' . In fact, if we again make the approximation that p « (p) to lowest order, then this residual 
baroclinic term is just proportional to (Vp'x VP'). Other contributions to Q include the nonlinear advection of the 
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mean meridional flow (the term involving (u m ), which is likely negligible as noted above), the Reynolds stress (the 
term involving u'), the Lorentz force (the term involving B), and the viscous diffusion (the term involving T>). 

We use the term turbulent stresses loosely, since it is clear that J- and Q may in principle include contributions from 
the large-scale Lorentz force and the viscous diffusion, which need not be turbulent. However, these are not likely to be 
important in the solar NSSL. The small molecular viscosity of the solar plasma makes viscous diffusion negligible and 
the persistence of the NSSL throughout the solar cycle suggests that it is not maintained by the large-scale Lorentz 
force. Thus, we expect that the dominant components of T and Q are indeed the turbulent Reynolds and Maxwell 
stresses associated with small-scale convection. 

APPENDIX B: AN ANALYTIC ILLUSTRATION OF GYROSCOPIC PUMPING 

Because it is essentially non-local and thus sensitive to the global geometry, boundary conditions, and inhomo- 
geneities, gyroscopic pumping is a subtle phenomenon that generally requires numerical calculations to find equilib- 
rium states. Still, analytic solutions can be f ound for ideal i zed ca s es. Here we present s uch a n analytic solution in 
the co ntext of the solar NSSL. For o t hers, see lHavnes et al.l (|1991f ). IGaraud fc Brummelll (2008), Garaud fc Arreguml 
IpOOl ). and IGaraud fc Bodenheimen (puTol ). 

We emphasize again the point made in i)2.f \ that the gyroscopic pumping equation @ can induce a meridional 
flow with relatively little impact on the rotational shear. Here we illustrate this by independently specifying both a 
differential rotation profile Q(X) and a zonal force T and then proceeding to derive a solution for the meridional flow 
that links the two. 

We restrict our attention to the northern hemisphere. Solutions for the southern hemisphere then follow by symmetry. 
Consider a simple, cylindrical angular velocity profile given by 

Q = Q a + AQ — . (Bl) 
R 

This, by construction, satisfies the steady state {d (w^) /dt = 0) meridional force balance (eq. ([8])) with B = Q = 0. In 
this appendix we will find a steady solution ^(A,z) to the zonal force balance equation ([2]) for the fi profile in (|BI[) 
and a specified torque J- . We will treat the differential rotation Af2 as a free parameter. Note that this is not a unique 
solution in the sense that other mean flow profiles 17(A) , ^(A, z) can be realized with the same J- ' . The equilibrium 
solution that will be realized in practice will depend on the initial and boundary conditions of the full time-dependent 
system. 

Here we consider a zonal torque of the form 

T =T a {r 2 +ar + b)r sin 2 9 cos 2 9 (Region I: r s < r < R) (B2) 

(Region 2: r < r s , z > z e ) (B3) 

F e {r,6) (Region 3: z < z e ) (B4) 

where Fq is the amplitude of the force and r s is the bottom of the NSSL (Region I). The coefficients a and b will be 
chosen to ensure that the meridional flow is continuous across r s . Also, we'll choose a, 6, and JFq to give us a negative 
J- in Region 1. The latitudinal dependence in (|B2I) is chosen to produce a vg that peaks at mid-latitudes, going as 
sin 9 cos 9 in the limit Aft —> (see eq. (|BI8|) below). The boundary layer at the equator (Region 3, equation (|B4[) ] is 
included to close the circulation cell in the northern hemisphere and to ensure that the net torque J v FdV — 0. Since 
is it a passive response, we will compute it only after we have computed the mean flow in Regions 1 and 2. The free 
parameters that specify T are thus the amplitude in the NSSL, J-q, and the locations of the boundary layers, r s and 

Z e - 

From equation (|BI[) we have 

/ r 

— = 2tt A(I + 7 A) (B5) 

where 7 = 3An/{2il R) . 
Equation yields the A component of the flow in Region 1: 

„ r 2 + ar + b Xz 2 _ . , . 

(pv x ) - /3 —— r (Region I) (B6) 



where (3 = J c q/(2Ho). In obtaining (|F36[) we have used the relations A = rsin# and z = rcos9. 
Equation © then yields 

*(\,z')= T l± x [I(\,z)-I(\,z b )] (z s <z'<z b ) (B7) 

*(A,z')= *(A,z s ) (z e <z'<z s ) 

*(A,z')= *s(A,z) (0<z'<z e ) 

where z b — \/ R 2 — A 2 as in Sj2j z s — \]r 2 s — A 2 , and we will specify ^3 (A, z) below. The solution in Region 1 involves 
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a 




Fig. 11. — The analytic solution for (a) C and (6) if according to equations ||B1[| and ||B7|I . Here the NSSL is a bit wide, r s = 0.90R 
(indicated by the dotted line), in order to highlight the structure of the boundary layer. Here we've taken 7 = 0.45-R - 1 and z e = 0.15-R. 
Contours in (ft) represent streamlines of the mass flux, with poleward flow at the surface. Dotted lines delineate Regions 1, 2, and 3. 

the following integral 
I(X,z)-- 



(r 2 +ar + b) — dz = — (r 2 + 2ar - 2b) + ( b 



A 2 
2 



In \z 



aX tan 



"(!) 



(B8) 



We have dropped the integration constant because we'll use this as a definite integral in what follows. 
In order to compute the coefficients a and b we will need the axial flow component, which follows from equation ([3]) 



(pvz) = -j8 



2 + 7A 
(I+7A) 1 



[I(X,z) 



d d 
-I(X,z)--I(X,z b ) 



(Region 1) 



(B9) 



with (pv z ) = (pv z ) \ z=Zs in Region 2. 

Note that the term in equation (|B9|) involving I(X, z) — I(X, z b ) vanishes at z — z b but as we'll see below, the term 
involving their derivatives does not. This is consistent with the impenetrable boundary condition v r = at r = R 
since the nonzero ve has both v\ a nd v z components. 

The derivatives in equation (IB9[) are given by 



d Xz 

I(X, z) = (r 2 + 2ar + 2b) —r - X In \z + r I 

2r 4 



2b -X 2 X 



and 



dX 



Tx I(X,z b ) 



z + r 2r 



— a tan 



(2b- R z 



2Rz b 



A In \z b + R\ 



2b -X 2 X 



a tan 1 I - 1 



Zb + R 2zb 



(BIO) 



(Bll) 



where we have used dzt/dX — —X/zj,. 

Now we proceed to compute the coefficients a and b. Since v\ = for r < r s , then continuity of v\ requires that it 
should vanish at r = r s as well, along with the forcing T . Using (|B2[) . this can be achieved if b — —r s (r s + a). 

As a consequence of the way we set up the problem, we also have 



Tx {Xpvx) = Yz {pVz) 







(Region 2) 



(B12) 



as well. It is 



So, in order for the circulation to be continuous, we require these expressions to hold at r = 
straightforward to show from equation (|B6j) that equation (|B12j) is satisfied at r = r s if r 2 + 2ar s + 36 = 0. Combining 
this with the previous expression then yields a = —2r s and b = r 2 . 

Note that with these values of a an d b, we have r 2 + ar + b = (r — r s ) . So, if we choose a negative value for Fq 
(which implies f3 < 0), then equation (|B2[) ensures that J- < in the NSSL, as desired. 
This completes the solution for \&(A, z) in Regions 1 and 2. For Region 3 we can then set 



Z (2z e -z) *(A,z e ) + — (z 



z e )*'(A,z e ) 



(B13) 



where vp' = /dz = (pv\). The value at z = z e is zero for r e < r 8 and is given by equation (|B6I) for r e > r s . Here 
r e = (z 2 + X 2 ) 1 / 2 . It is straightforward to show that this yields $ = at z = 0, implying no flow across the equatorial 
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plane, and furthermore, that d^>/dX, and d^/dz are continuous across z — z e . The analytic solution is plotted in 
Figure [TTJ 

The corresponding cylindrically outward flow in Region 3 is given by equation J5| 

(pv x ) = \{z e - z) *(A, z e ) + 2Z ~ Ze &(\, z e ) (Region 3) , (B14) 

z e z e 

while equation ([2]) gives the required torque 

Fz{\z) = ^{pv x ) (B15) 

with dC/d\ given by (|B5|) and (pv\) given by ()B14|) . Note that the volume- integrated torque is zero by construction, 
as a consequence of the impenetrable boundaries and equation ([2]) 

/ TdV = f ( P v m )-V£dV= f V-({pv m )£)-dS = . (B16) 
Jv Jv Js 

The axial flow in Region 3 follows from equations ([5]), (|B9[) . and (|B13|) 

(pv z ) = (2z e - z) (pv z ) \ z=Zs (Region 3, r < r s ) 

(pv z ) = (2z e - z) (pv z ) \ z=Ze - /J^^f {z - z e ) + f - - (Region 3, r > r s ) . 

It is instructive to verify that the flow at the outer boundary is indeed poleward. Substituting z = z\, and r = R 
into equations (|B6j) and (|B9j) and expressing the result in spherical coordinates yields 



BD 2 BD 2 
{PVX ) = TT ^-^ sin<W* , and (pv z ) = -^^-^ sin^cos* (r = R) (BIT) 

where D = R — r s is the thickness of the NSSL. Furthermore, 

i \ ■ n i \ n i \ FoD 2 sin cos , , ni „. 

(pvg) = -sme(pv z ) + cos9(pv x ) = —— ,. /Aa/n , . a (r = R) , BI8 

and (pv r ) = cos^ (pv z ) + sm9 (pv\) = at r = R. For < 0, equation (|Bl8|) implies a poleward flow. 

This is clearly a highly idealized depiction of a star. In the solar convective envelope, strong turbulent stresses T 
(and possibly Q) and baroclinic forcing B maintain a substantial non-cylindrical differential rotation profile (Fig. [l| 
and this will in turn determine how the meridional circulation streamlines close in Regions 2 and 3. In particular, 
the initial circul ations established by co nvection on a dynamical time scale are unlikely to penetrate much below the 
convection zone (|Gilman fc Mie sch 2004), although very weak gyroscopically-pumped circulations driven by convection 
will burrow downw ard on a radiative diffusion time scale, eventually reaching the deep radiativ e interior over the course 
of billions of years ( Spie gel fc Zahn|[l992HGough fc Mclntvrdll99 8; Gara ud fc Bru mmcll 2008). Here turbulent stresses 
are weak T rs and the rotation is nearly uniform so the circulation contours would follow cylindrical surfaces as in 
Region 2 of the present example. In any case, the solution discussed in this section is only intended to give a rough 
feel for how gyroscopic pumping in the NSSL might impact the dynamics in the upper convection zone. 

Still, equations (|Bl[) . (|B7[) . (|B8I) . and (|BI3j) . provide a steady, analytic solution to the continuity, momentum, and 
energy equations in the barotropic (or isentropic) limit (P — P{p)) under the influence of specified turbulent stresses 
analogous to the NSSL (J 7 < for r > 0.95). They also highlight the need for turbulent or baroclinic stresses in the 
meridional plane, B and/or Q, in order to maintain an axial rotational shear dfl/dz =/= 0. 

APPENDIX C: UNCURLING THE ZONAL VORTICITY EQUATION 

As described in $4.3[ we wish to derive the acceleration A m corresponding to the meridional stress Q. To do this, 
we must solve equation (|23[) for £. The expansion for O in equation (| L6[) implies that we can write £ as 



C(r,0) =sin0^Z„(r)cos n (n = l, 3, 5, T, 9). (CI) 



Each of the Z n (r) satisfy a recursive equation of the form 

K + \ Z 'n - ^ (« 2 + 3n + 2) Z n = -r„ - 1 (n + l) (n + 2) Z n+2 (C2) 

where primes denote (ordinary) derivatives with respect to r and Z n = for n even and for n > 9. The T n are given 
by 

ri = -2r^ e di r 3 = -2r (Sl e d 3 + Q. 2 di) T 5 = -2r {n e d 5 + fl 2 d 3 + Sl^di) (C3) 

T 7 = -2r {n 2 d 5 + n^ds) r 9 = -2r^ 4 d 5 (C4) 
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where 

2 2 4 2 4 

dx = n' e + -02 ^3 = ^2 - + -^4 d 5 = ft 4 + -fi 2 - -0 4 . (C5) 

The radial derivatives in equations (|C2|) and (IC5|) are discretized using a second-order finite difference scheme and the 
resulting matrix equations are solved subject to the boundary conditions Z n = at r = R and rZ' n + Z n = at r = r s . 
These boundary conditions ensure that the vertical acceleration vanishes at the surface (A r = at r = R) and the 
horizontal acceleration vanishes at the base of the NSSL (Ag = at r = r s ). 

APPENDIX D: DIFFUSIVE ANGULAR VELOCITY PROFILES 

In this Appendix we seek angular velocity profiles for which the net axial torque vanishes J- = 0, subject to the 
boundary conditions 

Q(r s ,6) = £l + fi 2 cos 2 + f2 4 cos 4 6» = uj + w 2 -P2(cos0) + LJiP^cosO) (Dl) 

and 

= (D2) 



on 

dr 



with coefficients 



at a specified matching layer, r = r s . Here S7o, f^2 and f^4 are fitting coefficients corresponding to the solar rotation 
profile at the mat chin g layer (r = r s ), as inferred from global helioseismology (Fig. [T]). The expression on the far 
right of equation (|Dlj) provides an alternate representation of f2(r s , 9) in terms of Legendre Polynomials P n (x). The 
Legendre coefficients are given by 

Cl Cl^ 2 4 8 

uj = Cl a + — + — uj 2 = 7:^2 + 3^4 w 4 = —^4 ■ (D3) 

3 5 3 7 35 

We choose r s to correspond to the location where the radial gradient of the spherically-averaged angular velocity 
profile passes through zero. This yields values of r s = 0.946-R, f2o/(27r) = 468 nHz, fl 2 /(2ir) = - 62.5 nHz, and 
A 4 /(2tt) = - 77.9 nHz. 

We consider two hypothetical forms for J 7 , representing turbulent diffusion and the mixing of angular momentum 
as discussed in <j5] In the first case the angular momentum flux is given by equation (|28|) . The assumption of no 
net torque and a constant density- weighted diffusion coefficient pv t then yields equation (f3~Tj): V(A 2 Vil„) = 0. The 
solution can be written as a Legendre series 

n v {r,8) = W a (r) + W 2 (r)P 2 (cos6) + W4(r)P 4 {cos6) (D4) 

W A {r) = 01 (r 4 + a 2 r" 7 ) (D5) 
W 2 (r) = b ir 4 + b 2 r- 7 + b 3 r 2 + b 4 r- 5 (D6) 
Wo(r) = cir 4 + c 2 r~ 7 + c 3 r 2 + c 4 r" 5 + c 5 r~ 3 + c 6 (D7) 

4 11 w 4 . 5 5 , . 

«2 = -zr m ax = -j— h = -ax b 2 = - ai a 2 (D8) 

5 9 2 

63 = y^ 2 r m 2 - -b ir 2 m + -b 2 r~ 9 b A = uj 2 r 5 m - b x r 9 m - b 2 r~ 2 - b 3 r 7 m (D9) 

Cl = ^4- C2 = ^4— C3 " C4 = ¥ (D10) 

c 5 = ^ (4c 1 r 7 „ - 7c 2 r m 4 + 2c 3 r£, - 5c 4 r m 2 ) c 6 = w - cir 4 ^ - c 2 r~ 7 - c 3 r 2 rn - c 4 r~ 5 - c 5 r" 3 . (Dll) 

For the second case we consider, that of angular momentum mixing, the angular momentum flux is given by equation 
(|2"9"|) and the rotation profile for T = and pv a = constant is given by equation (|3"2")l . V 2 £ a = 0. The solution is 
readily obtained through a Legendre series 

r n 1 a r — (n+1) 

£ B (r,0) = ^ P " Z7— T) L n P n (cos8) (n = 0, 2, 4, 6). (D12) 

n 1 m < Pn 1 m 

Here L n are the coefficients corresponding to the matching layer 

L 6 = - ^ for*, L 4 = ^ (0 4 - fo) r 2 „ + (D13) 



where 
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L 2 = \ (Aa - Clo) r 2 m + \U - f £ 6 L = il r 2 m + ^ - ~L 4 + ^L 6 (D14) 

and fi n — [(n — 2)/(n + 3)] r^ +1 . The corresponding rotation profile is then il a = \~ 2 £ a - 
The two solutions derived here, Sl v (r,8) and Vl a (r,9) are shown in Figure |H1 



